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Direct  numerical  simulations  of  the  interaction  of  a  premixed  flame  with  subsonic,  high-speed,  homogeneous,  isotropic,  Kolmogorov-type 
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unsplit,  reactive-flow  code.  A  simplified  reaction-diffusion  model  represents  a  stoichiometric  H2-air  mixture  under  the  assumption  of  the  Lewis 
number  Le  =  1.  Global  properties  and  the  internal  structure  of  the  flame  were  analyzed  in  an  earlier  paper,  which  showed  that  this  system 
represents  turbulent  combustion  in  the  thin  reaction  zone  regime  with  the  average  local  flame  speed  equal  to  its  laminar  value,  S j.  This  paper 
shows  that:  (1)  Flamelets  inside  the  flame  brush  have  a  complex  internal  structure,  in  which  the  isosurfaces  of  higher  fuel  mass  fractions  are 
folded  on  progressively  smaller  scales.  (2)  Global  properties  of  the  turbulent  flame  are  best  represented  by  the  structure  of  the  region  of  peak 
reaction  rate,  which  defines  the  flame  surface.  (3)  The  observed  increase  of  Sj  relative  to  Sj  exceeds  the  corresponding  increase  of  the  flame 
surface  area,  AT,  relative  to  the  surface  area  of  the  planar  laminar  flame,  on  average,  by  =  30%  and  occasionally  by  as  much  as  50%  in  the  course 
of  system  evolution.  This  exaggerrated  response  of  Sj  shows  that  Damkohler’s  paradigm  breaks  down  for  sufficiently  high-intensity  turbulence, 
namely  at  Karlovitz  numbers  Ka  i>  20,  even  in  the  flows  characterized  by  Le  =  1.  (4)  The  breakdown  is  the  result  of  tight  flame  packing  by 
turbulence,  which  causes  frequent  flame  collisions  and  formation  of  regions  of  high  flame  curvature  >1/5 j,  or  “cusps,”  where  SL  is  the  thermal 
width  of  the  laminar  flame.  (5)  The  local  flame  speed  in  the  cusps  substantially  exceeds  its  laminar  value,  which  results  in  a  disproportionately 
large  contribution  of  cusps  to  Sj  compared  with  the  flame  surface  area  in  them.  (6)  Results  suggest  the  existence  of  two  distinct  regimes  of  flame 
evolution:  the  linear  regime  present  at  lower  turbulent  intensities  when  Sj  Aj ,  and  the  nonlinear  regime,  which  is  dominated  by  the  contribu¬ 
tion  of  cusps  to  Sj  and  in  which  Sj  becomes  a  nonlinear  function  of  Aj-.  (7)  Finally,  the  criterion  for  transition  to  this  nonlinear  regime  is 
established,  and  key  properties  and  implications  of  this  regime  are  discussed. 
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The  Interaction  of  High-Speed  Turbulence  with  Flames:  Turbulent  Flame  Speed 


Abstract 

Direct  numerical  simulations  of  the  interaction  of  a  premixed  flame  with  subsonic,  high-speed,  homogeneous, 
isotropic,  Kolmogorov-type  turbulence  in  an  unconfined  system  show  anomalously  high  turbulent  flame  speeds,  St- 
Data  from  these  simulations  are  analyzed  to  identify  the  origin  of  this  anomaly.  The  simulations  were  performed 
with  Athena-RFX,  a  massively  parallel,  fully  compressible,  high-order,  dimensionally  unsplit,  reactive-flow  code. 
A  simplified  reaction-diffusion  model  represents  a  stoichiometric  H2-air  mixture  under  the  assumption  of  the  Lewis 
number  Le  =  1.  Global  properties  and  the  internal  structure  of  the  flame  where  analyzed  in  an  earlier  paper,  which 
showed  that  this  system  represents  turbulent  combustion  in  the  thin  reaction  zone  regime  with  the  average  local  flame 
speed  equal  to  its  laminar  value,  S  l-  This  paper  shows  that:  (1)  Flame  lets  inside  the  flame  brush  have  a  complex 
internal  structure,  in  which  the  isosurfaces  of  higher  fuel  mass  fractions  are  folded  on  progressively  smaller  scales. 
(2)  Global  properties  of  the  turbulent  flame  are  best  represented  by  the  structure  of  the  region  of  peak  reaction  rate, 
which  defines  the  flame  surface.  (3)  The  observed  increase  of  St  relative  to  Si  exceeds  the  corresponding  increase 
of  the  flame  surface  area,  At,  relative  to  the  surface  area  of  the  planar  laminar  flame,  on  average,  by  »  30%  and 
occasionally  by  as  much  as  50%  in  the  course  of  system  evolution.  This  exaggerrated  response  of  S  r  shows  that 
Damkohler’s  paradigm  breaks  down  for  sufficiently  high-intensity  turbulence,  namely  at  Karlovitz  numbers  Ka  >  20, 
even  in  the  flows  characterized  by  Le  =  1.  (4)  The  breakdown  is  the  result  of  tight  flame  packing  by  turbulence, 
which  causes  frequent  flame  collisions  and  formation  of  regions  of  high  flame  curvature  >  1  /6l,  or  “cusps,”  where 
8l  is  the  thermal  width  of  the  laminar  flame.  (5)  The  local  flame  speed  in  the  cusps  substantially  exceeds  its  laminar 
value,  which  results  in  a  disproportionately  large  contribution  of  cusps  to  St  compared  with  the  flame  surface  area  in 
them.  (6)  Results  suggest  the  existence  of  two  distinct  regimes  of  flame  evolution:  the  linear  regime  present  at  lower 
turbulent  intensities  when  St  At,  and  the  nonlinear  regime,  which  is  dominated  by  the  contribution  of  cusps  to  S  y 
and  in  which  S  t  becomes  a  nonlinear  function  of  /V/  .  (7)  Finally,  the  criterion  for  transition  to  this  nonlinear  regime 
is  established,  and  key  properties  and  implications  of  this  regime  are  discussed. 

Key  words:  Turbulent  premixed  combustion.  Turbulence,  Flamelet,  Turbulent  flame  speed.  Hydrogen 


1.  Introduction 

One  of  the  fundamental  questions  of  combustion  research  concerns  our  ability  to  understand  and  predict  the  rate 
of  energy  release,  or  equivalently  the  speed,  of  the  turbulent  flame,  which  is  typically  the  main  factor  controlling  the 
evolution  and  dynamics  of  the  overall  system.  To  achieve  this,  two  key  aspects  of  the  combustion  process  must  be 
determined:  (1)  its  local  properties,  namely  the  local  speed  of  flame  propagation,  and  (2)  its  global  characteristics, 
i.e.,  the  overall  structure  of  the  turbulent  flame  that  connects  the  local  burning  velocity  at  each  point  of  the  flame 
surface  with  the  turbulent  flame  speed.  In  general,  such  local  and  global  characteristics  are  not  universal  for  the 
combustion  process.  They  can  vary  substantially  both  in  space  and  time  due  to  the  unsteadiness,  inhomogeneity,  and 
anisotropy  of  the  turbulent  flow  associated  with  the  system  geometry,  presence  of  walls  and  boundaries,  change  in  the 
flow  conditions,  etc. 

A  remarkably  simple  yet  powerful  paradigm,  which  suggests  a  description  both  of  the  local  flame  properties  and 
of  their  connection  with  the  turbulent  flame  speed,  was  proposed  by  Damkohler  almost  70  years  ago  [1]  (also  see 
reviews  by  [2,  3]).  Today  it  still  remains  a  cornerstone  of  our  understanding  of  turbulent  combustion.  This  paradigm 
is  formulated  for  a  turbulent  flame  that  is,  on  average,  a  planar,  quasi-one-dimensional  structure  propagating  in  a 
steady,  homogeneous,  and  isotropic  turbulent  flow. 

According  to  the  Damkohler’s  suggestion,  the  local  flame  speed,  S  /,  is  determined  only  by  the  intensity  of  turbulent 
motions  on  scales  smaller  than  the  laminar  flame  width,  6l-  If  that  intensity  is  sufficiently  low,  the  flame  propagates 
at  each  point  of  its  surface  with  the  laminar  flame  speed,  S  l-  On  the  other  hand,  if  small-scale  motions  are  energetic 
enough,  they  are  able  to  penetrate  the  internal  structure  of  the  flame  and  broaden  it.  Thus  molecular  diffusivity 
becomes  enhanced  by  the  turbulent  transport  associated  with  such  small-scale  turbulence,  thereby  increasing  the  local 
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flame  speed,  so  that  S /  ~  5 /.(/>///?) I,/2,  where  D  and  />/  are,  respectively,  the  molecular  and  the  effective  turbulent 
diffusivities. 

At  the  same  time,  the  overall  structure  of  the  turbulent  flame  is  determined  by  large-scale  turbulent  motions 
which  stretch  and  fold  the  flame,  thereby  increasing  its  surface  area.  At.  Consequently,  since  at  each  point  the  flame 
propagates  with  the  speed  5  /,  the  increase  of  the  turbulent  flame  speed,  S  r,  relative  to  that  of  the  planar  laminar  flame 
is  equal  to  the  increase  of  the  flame  surface  area,  i.e.,  S -/  /.S' /  =  At/At,  where  A/  is  the  surface  are  of  the  planar 
laminar  flame.  Therefore,  in  this  paradigm,  in  order  to  determine  S l,  it  is  sufficient  to  model  only  one  aspect  of  the 
process  of  flame  wrinkling  and  folding  by  turbulence,  namely  the  resulting  increase  in  the  flame  surface  area. 

This  is  essentially  a  kinematic  model  which  neglects  any  potential  feedback  of  the  flame  on  the  turbulent  flow.  At 
its  root  is  one  key  assumption:  S  /  is  determined  only  by  scales  A  <  6l,  and  it  is  constant  throughout  the  flame  surface. 
This  implies  that  S  /  cannot  depend  on  the  local  flame  geometry,  i.e.,  flame  curvature,  which  varies  between  different 
points  of  the  flame  surface.  Consequently,  S  /  cannot  depend  on  scales  A  >  6/  that  determine  such  local  geometry. 
Damkohler’s  concept  is  thus,  in  effect,  a  local  model,  since  in  it  the  long-range  velocity  correlations  do  not  affect 
the  local  properties  of  the  combustion  process.  This  assumption  then  allows  one  to  reduce  the  full  complexity  of  the 
three-dimensional  (3D)  configuration  of  the  turbulent  flame  to  one  very  simply  measure  -  flame  surface  area.  The 
result  is  a  remarkably  simple  way  to  connect  the  local  flame  properties  with  the  global  characteristics  of  the  turbulent 
flame,  namely  its  speed. 

The  limitations  of  this  assumption  were  known  practically  at  the  time  when  this  model  was  proposed  [4].  When 
the  Lewis  number  Le  +  1,  the  local  flame  speed  can  vary  with  the  curvature  of  the  flame.  In  particular,  if  thermal 
conduction  is  stronger  than  molecular  diffusion,  S  /  will  increase  in  concave  regions  of  the  flame  and  decrease  in 
convex  regions.  As  a  result,  scales  A  >  6l  also  affect  5  /  by  stretching  and  folding  the  flame.  These  effects  of  flame 
strain  and  curvature  can  then  be  incorporated  into  the  relation  between  At  and  S  by  means  of  the  stretch  factor  I 
([3],  also  see  [5]  for  the  review  of  the  theory  of  flame  stretch), 
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In  principle,  at  this  point  the  simplicity  of  the  description  of  the  turbulent  flame  structure  is  lost,  since  I  requires 
knowledge  not  just  of  the  flame  surface  area  but  of  how  the  flame  is  actually  folded  inside  the  flame  brush.  Neverthe¬ 
less,  for  the  Le  —  1  reactive  mixtures,  /  by  definition  must  be  unity,  and  the  details  of  the  flame  configuration  should 
not  affect  the  local  speed  of  flame  propagation. 

Despite  the  broad  acceptance  of  this  picture,  there  are  only  a  few  studies  aimed  at  directly  verifying  it.  From  an 
experimental  point  of  view,  simultaneously  determining  At  and  .S'  y  in  a  highly  turbulent  flow  with  sufficient  accuracy 
is  associated  with  substantial  practical  difficulties.  Resulting  uncertainties  typically  do  not  allow  robust  verification  of 
the  balance  between  these  two  quantities. 

Several  numerical  studies,  however,  do  consider  this  issue.  Bell  et  al.  [6]  analyzed  this  balance  in  the  3D  direct 
numerical  simulations  (DNS)  of  the  statistically  planar  turbulent  methane  flame  and  determined  that  /  deviates  from 
unity  by  <  10%.  A  similar  estimate  was  obtained  by  Khokhlov  [7]  in  the  3D  simulations  of  the  Rayleigh-Taylor 
driven  thermonuclear  flames  in  degenerate  matter.  Both  of  these  studies,  however,  considered  reactive  mixtures  in 
which  Le  +  1  and,  therefore,  /  would  not  be  expected  to  be  exactly  equal  to  one.  In  this  context,  the  result  of  [7]  is 
particularly  interesting  since  degenerate  matter  is  characterized  by  Le  »  1  and  typically  in  excess  of  103.  Both  studies, 
however,  did  consider  fairly  low  turbulent  intensities,  which  precluded  the  formation  of  complex,  tightly  folded  flame 
configurations.  For  instance,  in  [6],  the  characteristic  turbulent  velocity  was  «4.3Sl,  while  in  [7],  it  reached  values 
<  1 2, S'  /, .  Hawkes  &  Chen  [8]  analyzed  the  balance  between  At  and  .S'  y  in  the  simulations  of  lean  statistically  flat 
turbulent  methane  flames  at  higher  turbulent  intensities  <  28.55 l-  They  found  the  values  of  /  generally  within  a  few 
percent  of  unity,  with  the  exception  of  the  methane-hydrogen  flame  for  which  I  =  1.18  was  determined.  Interpretation 
of  the  results  of  [8],  however,  is  complicated  by  the  two-dimensional  nature  of  their  simulations  and  by  the  fact  that 
the  turbulent  integral  scale  was  <  d/,,  which  substantially  suppressed  flame  wrinkling  and  prevented  the  formation  of 
the  tightly  packed  flame  structure. 

These  studies  generally  support  Damkohler’s  concept  by  showing  that  the  increase  in  5  r  is  heavily  dominated  by 
the  growth  of  At  with  the  effects  of  the  flame  configuration  inside  the  flame  brush  playing  a  minor  role.  They  do  not, 
however,  answer  the  following  two  questions.  First,  while  I  appears  to  be  close  to  unity  for  the  Le  +  1  mixtures,  is 
/  indeed  exactly  equal  to  one  in  the  case  of  Le  —  1?  Since  realistic  reactive  mixtures  typically  do  not  have  Le  -  1 
exactly,  this  question  should  be  understood  in  the  sense  whether  /  — >  1  as  Le  — >  1.  Second,  how  does  the  value  of 
/  change  with  turbulent  intensity?  In  other  words,  does  there  exist  a  regime  in  which  Damkohler’s  concept  can  be 
expected  to  break  down,  leading  to  values  of  /  substantially  different  from  unity  not  only  for  Le  +  1,  but  also  for 
Le  -  1  systems? 

The  last  question  is  of  particular  importance  since  the  existence  of  such  a  regime  would  also  imply  the  breakdown 
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of  flamelet  models  that  focus  on  determining  A/  or  its  volumetric  equivalent,  flame  surface  density  Zy,  in  order  to 
predict  S  r  ■  Furthermore,  in  such  a  regime  the  local  flame  speed  would  no  longer  be  determined  only  by  the  small 
scales  A  <  6l  but,  rather,  by  the  full  turbulent  spectrum.  This  would  mark  the  breakdown  of  the  locality  inherent  in 
Damkohler’s  model.  In  particular,  flame  propagation  in  this  regime  could  no  longer  be  considered  as  a  local  process 
on  any  scale,  since  the  local  quantity  S  /  would  depend  on  the  nonlocal  long-range  velocity  correlations  that  potentially 
span  the  full  size  of  the  system. 

The  objective  of  this  paper  is  to  begin  addressing  these  two  questions  by  considering  the  balance  between  Aj  and 
S  t  in  the  ideally  infinite,  statistically  planar,  turbulent  flame  interacting  with  the  high-speed,  steady,  homogeneous, 
isotropic  turbulence.  We  focus  on  the  Le  =  1  situation  to  exclude  thermodiffusive  effects  as  a  potential  source  of  the 
non-unity  values  of  I. 

This  paper  continues  the  analysis  of  the  simulations  first  presented  by  us  in  [9],  These  calculations  model  flame 
interaction  with  turbulence  of  sufficiently  high  intensity  to  represent  the  regime  that  is  borderline  between  thin  and 
broken  reaction  zones,  according  to  the  traditional  combustion  regime  diagrams  [10,  11,  12]  (cf.  Fig.  2  in  [9]).  In  other 
words,  we  consider  the  fastest  turbulence  which  has  been  hypothesized  to  allow  the  existence  of  flamelets  with  the 
internal  structure  of  the  reaction  zone  essentially  unaffected  by  the  turbulent  transport.  The  turbulent  r.m.s.  velocity  of 
the  flow  in  cold  fuel  is  Urms  ~  3 5 .S' / ,  leading  to  the  Damkohler  number  Da  =  0.05  and  Gibson  scale  Lg  ~  3  x  1 0  46l 
(Table  2). 

The  primary  focus  of  [9]  was  a  detailed  study  of  flame  properties  and  evolution  in  the  presence  of  such  high-speed 
turbulence.  In  particular,  a  method  was  presented  which  allowed  the  direct  determination  of  the  internal  structure 
of  the  flame  based  on  its  actual  3D  configuration  inside  the  flame  brush.  The  analysis  showed  that  the  preheat  zone 
is  broadened  by  turbulence  while  the  reaction-zone  structure  remains  virtually  identical  to  that  of  the  planar  laminar 
flame.  This  study  demonstrated  that,  even  under  the  action  of  such  intense  turbulence,  the  system  indeed  can  be 
robustly  classified  to  represent  the  thin  reaction  zone  regime. 

The  accurate  determination  of  the  reaction-zone  structure  in  [9]  showed  that  turbulence  is  not  able  to  penetrate  the 
flame  and  thus  does  not  enhance  diffusive  processes.  Consequently,  the  flame  must  propagate  locally  with  the  laminar 
speed  S / .  Given  that  we  are  considering  the  Le  -  1  system,  this  would  imply  that  St/S l  =  At /Al  with  a  high  degree 
of  accuracy,  according  to  Damkohler’s  concept. 

Before  the  relation  between  5  t  and  ,4y  can  be  studied,  however,  the  question  of  the  definition  of  the  flame  surface 
area  must  first  be  revisited.  This  question,  which  has  a  seemingly  simple  answer  at  low  turbulent  intensities,  becomes 
quite  non-trivial  in  the  case  of  a  high-speed  flow  such  as  the  one  considered  here.  The  flamelet  structure  determined 
in  [9]  is  the  manifestation  of  the  fact  that  the  ability  of  the  turbulent  cascade  to  penetrate  the  flamelet  interior  changes 
substantially  from  the  colder  parts  of  the  flamelet  in  the  outer  preheat  zone  to  the  hotter  regions  in  the  reaction  zone. 
One  prominent  consequence  of  this,  observed  in  [9],  was  wrinkling  of  the  turbulent  flame  on  much  smaller  scales  on 
the  fuel  side  than  on  the  product  side.  This  shows  that  effects  of  turbulent  motions  on  scales  A  <  6l  vary  significantly 
throughout  the  flamelet  interior,  causing  varying  response  of  its  different  parts  to  the  stretching  and  folding  action 
of  turbulence.  Therefore,  the  flamelet,  packed  into  the  flame  brush  by  intense  turbulence,  can  be  expected  to  have  a 
complex  internal  structure  with  the  isosurfaces  of  different  values  of  the  fuel  mass  fraction,  Y ,  or  temperature,  T ,  not 
being  parallel  to  each  other  and,  most  importantly,  not  having  the  same  surface  area. 

Our  starting  point  in  this  work  is  thus  the  following  question:  What  is  the  flame  surface  area  in  high-speed  turbulent 
flows?  Answering  this  will  then  allow  us  to  determine  whether  there  is  indeed  a  perfect  balance  between  the  increase 
of  the  burning  speed  of  the  turbulent  flame  and  its  surface  area,  as  suggested  by  Damkohler’s  concept,  when  Le  =  1. 

2.  Numerical  method  and  simulations  performed 

2.1.  Physical  and  numerical  models 

Here  we  summarize  the  physical  model,  the  numerical  method  used,  and  key  aspects  of  the  simulation  setup.  A 
detailed  discussion  can  be  found  in  [9], 

We  solve  the  system  of  unsteady,  compressible,  reactive  flow  equations. 


f  +  V(pu) 

=  0, 

(2) 

+  V  •  (pu  <8>  u)  +  VP 
at 

=  0, 

(3) 

^  +  V  •  ((£  +  B)u)  +  v  •  (/fvr) 

=  pqY, 

(4) 

^  +  V  •  (pTu)  +  V  •  (y oDVY ) 

=  PY. 

(5) 
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Here  p  is  the  mass  density,  u  is  the  velocity,  E  is  the  energy  density,  P  is  the  pressure,  Y  is  the  mass  fraction  of  the 
reactants,  q  is  the  chemical  energy  release,  and  Y  is  the  reaction  source  term.  The  coefficients  of  thermal  conduction, 
K,  and  molecular  diffusion,  D,  are 

rj-'n  j^n 

D  =  D0  — ,  — -  =  kq  — ,  (6) 

P  pCP  p 

where  Do,  kq,  and  n  are  constants,  and  Cp  =  -yR/M(y  -  1)  is  the  specific  heat  at  constant  pressure.  The  equation  of 
state  is  that  of  an  ideal  gas.  Chemical  reactions  are  modeled  using  the  first-order  Arrhenius  kinetics 

dY  (  Q\ 

-sY  =  -pYB^{-±\  (7) 

where  B  is  the  pre-exponential  factor,  Q  is  the  activation  energy,  and  R  is  the  universal  gas  constant. 

Table  1  summarizes  the  parameters  of  the  physical  model  used  as  well  as  the  resulting  properties  of  the  planar 
laminar  flame  [9],  These  parameters  are  based  on  a  simplified  reaction-diffusion  model  designed  to  represent  the 
stoichiometric  H2-air  mixture  [13]. 

Eqs.  (2)-(7)  are  solved  using  the  code  Athena-RFX  [9]  -  the  reactive  flow  extension  of  the  magnetohydrodynamic 
code  Athena  [14,  15].  Athena-RFX  is  a  fixed-grid,  massively  parallel,  finite-volume,  fully  conservative  code,  ft 
implements  a  variant  [15]  of  the  fully  unsplit  corner-transport  upwind  (CTU)  algorithm  [16]  and  its  3D  extension  pre¬ 
sented  in  [  1 7] ,  in  conjunction  with  the  PPM  spatial  reconstruction  [18]  and  the  approximate  nonlinear  HLLC  Riemann 
solver.  Overall,  the  code  achieves  3'^-order  accuracy  in  space  and  2"rf-order  accuracy  in  time.  A  detailed  description 
and  an  extensive  suite  of  tests  of  the  hydrodynamic  integration  algorithm  can  be  found  in  [14,  15].  Implementation 
of  the  reactive-diffusive  extensions  in  Athena-RFX,  along  with  the  results  of  tests  including  convergence  studies,  is 
discussed  in  detail  in  [9,  19]. 

The  simulations  presented  here  model  flame  interaction  with  steady,  homogeneous,  isotropic  turbulence,  described 
by  the  classical  Kolmogorov  theory  [20],  Turbulence  driving  is  implemented  using  a  spectral  method  [9,  19]  similar 
to  the  one  used  in  [21,  22],  In  this  method,  velocity  perturbations  du(k)  are  initialized  in  Fourier  space  with  each 
component  6u'( k)  being  an  independent  realization  of  a  Gaussian  random  field  superimposed  with  a  desired  energy 
injection  spectrum.  Subsequently,  non-solenoidal  components  of  du(k)  are  removed  to  ensure  that  V  •  c)u(x)  =  0. 
An  inverse  Fourier  transform  of  <5u(k)  gives  dii(x),  the  velocity  perturbation  field  in  the  physical  space.  The  c)ii(x)  is 
normalized  to  provide  the  constant  rate  s  of  kinetic-energy  injection,  and  the  total  momentum  in  the  perturbation  field 
is  subtracted  from  du(x)  to  ensure  that  no  net  momentum  is  added  to  the  domain,  i.e.,  f  p8u  =  0.  Resulting  velocity 
perturbations  are  added  to  the  flow  field  u(x)  on  every  time  step,  and  the  overall  perturbation  pattern  is  regenerated  at 
every  time  interval  A tvp  ~  5  A x/cs,  where  cs  is  the  maximum  sound  speed  in  the  domain. 

Energy  is  injected  only  at  the  scale  of  the  domain  width,  L,  to  obtain  the  Kolmogorov-type  spectrum.  The  resulting 
turbulence  is  statistically  steady,  isotropic,  and  homogeneous  with  the  inertial  range  of  the  energy  cascade  extending 
all  the  way  to  the  energy  injection  scale  (cf.  Fig.  1  in  [9]).  Moreover,  since  the  velocity  perturbation  field  is  divergence- 
free,  no  compressions  or  rarefactions  are  artificially  induced  as  a  result  of  driving. 

It  can  be  seen  in  eqs.  (2)-(5)  that  we  do  not  explicitly  include  physical  viscosity,  but  instead  we  rely  on  numerical 
viscosity  to  provide  the  kinetic-energy  dissipation.  By  systematically  varying  the  resolution  in  this  approach,  the 
energy  spectrum  can  be  extended  to  progressively  smaller  scales  A  «  5l  while  maintaining  the  spectrum  constant  on 
larger  scales.  This  allows  us  to  vary  the  intensity  only  of  the  small-scale  turbulent  motions,  and  thus  to  investigate 
their  effects  by  differentiating  them  from  the  effects  of  large  scales.  Such  analysis  also  shows  the  range  of  scales 
which  must  be  resolved  in  a  numerical  simulation  in  order  to  capture  accurately  the  evolution  of  the  turbulent  flame. 
This  is  particularly  important  in  the  context  of  high-speed  turbulent  reactive  flows,  in  which  it  is  typically  impossible 
to  resolve  the  Kolmogorov  scale.  The  analysis  presented  in  [9]  demonstrated  that  scales  A  «  6i  have  virtually  no 
effect  on  the  properties  of  the  turbulent  flame  with  the  exception  of  the  degree  of  flame  surface  wrinkling  on  the  fuel 
side.  This  showed  that  the  evolution  of  the  turbulent  flame  can  be  accurately  reproduced  without  the  need  to  resolve 
scales  A  <sc  6P.  We  refer  to  [9]  for  the  detailed  discussion  of  this  issue  and,  in  particular,  for  the  discussion  of  the 
applicability  of  the  obtained  results  to  the  actual  stoichiometric  Hj-air  mixture. 

2.2.  Summary  of  simulations 

Key  parameters  of  the  simulations  discussed  here  are  summarized  in  Table  2.  The  main  difference  between  the 
three  calculations  is  the  resolution,  which  progressively  increases  from  Ax  =  8l / 8  in  SI  to  Ax  =  6l/ 32  in  S3.  The 
domain  in  all  cases  was  initialized  with  uniform  density  po  and  temperature  To  (see  Table  1).  In  SI  and  S2,  initial  fluid 
velocities  were  set  to  0.  In  contrast,  the  velocity  field  in  S3  was  initialized  with  the  ideal  energy  spectrum  S(k)  oc  k  5/5 
extending  from  the  energy  injection  scale  L  to  the  numerical  Kolmogorov  scale  i]  =  2Ax.  This  initial  spectrum  was 
normalized  to  ensure  that  at  t  —  0  the  total  kinetic  energy  in  the  domain  was  equal  to  its  predicted  steady-state  value 
as  described  in  [19].  In  SI  and  S2,  the  flow  field  was  allowed  to  evolve  for  the  time  tjgn  =  3t(y/,  and  in  S3  for  the 
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time  tign  =  2t€Ci,  to  develop  the  steady-state  turbulent  flow  field.  At  the  time  t,x„ ,  a  planar  laminar  flame  with  its  front 
parallel  to  the  x—y  plane1  was  initialized  in  the  domain  with  the  values  of  p,  T ,  and  Y  based  on  the  exact  laminar  flame 
solution.  The  velocity  field  was  not  modified,  thereby  preserving  the  structure  that  developed  during  the  equilibration 
stage. 

Prior  to  tign,  all  domain  boundaries  are  periodic.  At  tlx„ ,  boundary  conditions  (BCs)  along  the  left  and  right  z- 
boundaries  are  switched  to  zero-order  extrapolation  in  order  to  prevent  pressure  build-up  inside  the  domain  and  to 
accommodate  the  resulting  global  fluid  flow.  Reflections  of  the  acoustic  perturbations  caused  by  such  BCs  can  be 
neglected  in  comparison  with  the  acoustic  noise  generated  by  the  turbulent  flow  field  itself,  and  overall  such  choice  of 
the  BCs  was  found  not  to  introduce  any  unphysical  artifacts  and  thus  not  to  affect  the  solution  accuracy.  Furthermore, 
any  potential  effect  of  the  outflow  BCs  was  minimized  by  the  large  length-to-width  ratio  of  the  computational  domain, 
which  allowed  us  to  keep  the  flame  brush  sufficiently  far  from  the  boundaries  at  all  times  [9].  Hereafter,  for  simplicity 
we  refer  to  the  moment  of  ignition  as  t  —  0. 

Energy  is  injected  into  the  domain  with  the  constant  rate  e  per  unit  volume  for  the  total  duration  of  the  simulations 
to  provide  steady  driving  of  the  turbulent  flow.  The  value  of  s  was  chosen  to  produce  a  high-intensity  turbulence  field 
that,  however,  was  weak  enough  to  minimize  the  probability  of  creating  weak  transonic  shocklets  arising  from  the 
intermittency  in  turbulent  flow.  The  characteristic  turbulent  velocities  and  the  turbulent  integral  scale  of  the  steady- 
state  flow  prior  to  the  moment  of  ignition  are  listed  in  Table  2.  Corresponding  values  of  the  Damkohler  number  and 
the  Gibson  scale,  given  in  the  table,  show  that  Da  <K  1  and  Lq  «  by ■  Prior  analysis  of  the  resulting  structure  of  the 
turbulent  flame  demonstrated  that,  despite  the  high  turbulent  intensity,  small-scale  turbulence  fails  to  penetrate  the 
internal  structure  of  the  flame  and  the  system  evolves  in  the  thin  reaction  zone  regime  [9], 

Detailed  discussion  of  the  turbulent  flame  evolution  in  the  simulations  is  given  in  [9] .  Here,  for  completeness,  we 
reproduce  in  Fig.  1  the  time  histories  of  the  flame-brush  width,  by,  normalized  by  by ,  as  well  as  Sr,  normalized  by 
S l,  (see  Fig.  4  in  [9])  for  all  three  calculations,  since  these  two  quantities  are  used  extensively  in  this  work. 

The  width  of  the  turbulent  flame  brush  is  defined  as 


5r 


—  -  I  .max  zq  jnin, 


(8) 


where  zoim;„  and  z\xnax  are  defined  as 


20 ,min  =  max(z)  :  Y(x,y,z)  <  0.05  V  (. x,y,z  <  z0,mm), 
Z\,max  =  min(z)  :  Y(x,y,z )  >  0.95  V  (x,y,z  >  zhmax). 


In  other  words,  z<;imin  marks  the  rightmost  .v-y-plane  to  the  left  of  which  is  pure  product,  while  z\  max  marks  the  leftmost 
.r-y- plane  to  the  right  of  which  is  pure  fuel.  This  is  illustrated  in  Fig.  2.  The  turbulent  flame  speed  is  defined  as 


Sr  - 


PoL1  ’ 


(10) 


where  Mr  is  the  total  rate  of  fuel  consumption  inside  the  flame  brush,  i.e.,  the  total  mass  of  reactants  which  is  trans¬ 
formed  to  product  per  unit  time.  The  detailed  discussion  of  this  choice  of  the  definition  of  S  y  can  be  found  in  [9]. 

Finally,  Table  3  lists  the  time-averaged  values  of  both  by  /  by  and  St/S  l  [9]  as  well  as  of  other  key  characteristics 
of  the  turbulent  flame  discussed  below.  Table  3  also  shows  the  corresponding  order  of  self-convergence  for  each 
quantity  listed.  Since  the  computational  cell  size  decreases  progressively  by  a  factor  of  2  for  each  simulation,  the 
order  of  self-convergence  0((p)  of  a  variable  (f>  is  defined  as 


O(f)  =  log2 


l<Ast  ~  <As3l 
\fsi  ~  0S3l, 


(ID 


3.  Flame  surface  area  and  its  relation  to  the  turbulent  flame  speed 

3.1 .  Surface  area  of  the  fuel  mass  fraction  isosurfaces 

Fig.  3a, b  shows  the  evolution  of  the  surface  areas  Ao.ot  and  A0.99  of  the  fuel  mass  fraction  isosurfaces  ot  Y  -  0.01 
and  Y  =  0.99,  respectively.  We  use  the  notation  Ay  =  A(Y  =  Y'),  and  all  isosurfaces  are  constructed  using  the 
“marching  cubes”  algorithm.  These  isosurfaces  represent  two  outer  boundaries  of  the  flame  brush,  separating  it  from 
pure  product  and  pure  fuel.  All  values  are  normalized  by  the  domain  cross-section  Lr  corresponding  to  the  surface  area 


1  The  longest  dimension  of  the  domain  is  along  the  2-axis. 
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of  the  planar  laminar  flame  unperturbed  by  turbulence.  In  principle,  the  normalization  should  be  performed  over  some 
average  surface  area  of  the  flame  brush,  thus  reflecting  its  overall  large-scale  shape.  With  sufficiently  good  accuracy, 
however,  the  turbulent  flame  in  the  system  considered  here  can  be  viewed,  on  average,  as  a  planar  propagating  front, 
and  thus  we  normalize  over  its  area  L2. 

Fig.  3a  shows  that  Ao.oi  exhibits  substantial  variability,  similar  to  c5y  and  St  (cf.  Fig.  1).  There  is  a  clear  corre¬ 
spondence  between  the  peaks  and  troughs  in  the  values  of  Ao.oi  and  S  y  and,  to  a  lesser  extent,  by.  Moreover,  in  a 
manner  similar  to  by  and  by,  /to  ol  shows  less  variability  with  increasing  resolution,  and  it  appears  on  average  to  have 
converged. 

Such  manifest  correlation  in  the  behavior  of  Ao.oi  and  by  is  the  reflection  of  the  evolutionary  cycle  of  the  flame 
brush,  discussed  in  [9],  It  was  shown  that  in  the  turbulent  flame,  which  on  average  is  in  a  steady  state,  the  influx 
of  fresh  fuel  and  the  rate  of  its  consumption  never  perfectly  balance  each  other.  Instead,  periodically  either  one  or 
the  other  process  dominates.  In  particular,  an  increase  in  the  flame  surface  area  inside  the  brush,  and  the  associated 
increase  in  by  ,  leads  to  a  higher  b  y .  This  causes  rapid  consumption  of  fuel  inside  the  brush,  decreasing  both  the  flame 
surface  area  and  the  overall  width  of  the  flame  brush.  As  a  result,  a  slower,  less  convolved  flame  emerges,  which  is 
thus  more  susceptible  to  the  action  of  turbulence.  This  leads  to  increased  folding  and  stretching  of  the  flame,  and  the 
cycle  repeats. 

The  behavior  of  A0.99  is  pronouncedly  different.  There  appears  to  be  no  correlation  between  variations  in  A 0.99 
and  either  6j  or  by.  Moreover,  with  increasing  resolution,  the  magnitude  of  the  variations  increases,  and  there  is  no 
evidence  of  convergence  of  the  growing  average  values  of  A0.99.  We  will  consider  the  correlation  between  A(Y)  and 
both  b t  and  by  in  a  more  quantitative  form  in  §  3.3. 

These  conclusions  regarding  the  change  of  Ao.oi  and  A0.99  with  resolution  are  supported  by  Fig.  3c,  which  shows 
the  full  normalized  time-averaged  distributions  A(Y)/L2  for  all  values  of  Y .  In  simulation  S 1 ,  Ao.oi  and  A0.99  are  nearly 
equal.  With  increasing  resolution,  however,  they  diverge  to  the  point  that,  in  S3,  there  is  almost  a  factor  of  2  difference 
between  them.  Values  of  Ao.oi  steadily  decrease  and  they  indeed  demonstrate  convergence.  At  the  same  time,  A0.99 
increases  with  resolution,  and  the  difference  between  S2  and  S3  is  only  marginally  smaller  than  between  S 1  and  S2. 

This  behavior  is  part  of  a  broader  qualitative  change  in  the  overall  shape  of  the  distribution  of  A  which  occurs 
around  Y  ~  0.5,  i.e.,  at  the  boundary  of  the  reaction  and  preheat  zones.  The  lower  resolution  of  SI  suppresses  small- 
scale  turbulent  motions,  which,  in  turn,  causes  less  flame  wrinkling  in  the  preheat  zone.  As  a  result,  the  flame-brush 
surface  appears  similar  both  on  the  product  and  fuel  sides  (cf.  Fig.  3  in  [9]).  With  higher  resolution,  however,  turbulent 
flame  wrinkling  becomes  more  pronounced  with  increasing  distance  from  the  reaction  zone,  and  the  distribution  of 
A  develops  a  distinctive  inverted-S  shape.  The  consequence  of  this  was  observed  in  [9]  (see  Fig.  3  therein),  which 
showed  a  much  more  highly  convolved  flame-brush  surface  on  the  fuel  side  in  calculations  S2  and  S3.  We  will  discuss 
the  role  of  small-scale  turbulence  in  further  detail  in  §  4.2. 

Profiles  of  A  in  S2  and  S3  are  closer  to  each  other  in  the  reaction  zone  than  in  the  preheat  zone.  Table  3  shows  that 
in  the  region  of  peak  reaction  rate,  i.e.,  at  Y  w  0.15,  A(Y)/L 2  exhibits  the  3,Y/-order  convergence.  This  is  substantially 
faster  than  what  would  be  expected  for  an  overall  2"d-order-accurate  code.  Moreover,  in  all  three  calculations,  A  varies 
the  slowest  inside  the  reaction  zone,  i.e.,  for  Y  «  0.15  -  0.5.  This  shows  that  isosurfaces  of  the  fuel  mass  fraction, 
on  average,  follow  each  other  closely  within  the  reaction  zone,  which  is  thus  folded  and  stretched  by  turbulence  as  a 
coherent  structure.  This  is  in  agreement  with  the  very  low  variability  in  the  reaction  zone  of  the  instantaneous  profiles 
of  Y  and  7’,  which  represent  the  internal  flamelet  structure  [9]. 


3.2.  Relation  between  S  y  and  A(Y) 


Results  presented  in  §  3.1  demonstrate  that  the  flamelet,  folded  inside  the  turbulent  flame  by  high- intensity  tur¬ 
bulence,  has  a  complex  internal  structure  with  a  distinctly  different  response  of  its  various  parts  to  the  action  of 
turbulence.  Consider  the  modified  eq.  (1) 


St  A(Y) 
Sl  O  ’ 


(12) 


in  which  the  surface  area  of  the  turbulent  flame  Ay  is  now  represented  by  the  isosurface  area  A(T),  and  we  set  the 
surface  area  of  the  planar  laminar  flame  Ay  equal  to  the  domain  cross-section  area  L2.  Using  data  for  A(Y)  and 
St,  eq.  (12)  then  allows  us  to  assess  the  extent  to  which  Damkohler’s  concept  is  valid  in  the  regime  modeled  in  our 
simulations.  Since  A  is  not  a  constant,  however,  but  rather  it  depends  on  the  choice  of  Y ,  I  now  also,  formally,  becomes 
a  function  of  Y. 

As  an  example,  Fig.  4a  shows  the  evolution  of  the  stretch  factor  / 0.15  =  I(Y  =  0.15)  based  on  the  area  of  the 
isosurface  that  corresponds  to  the  peak  reaction  rate  in  the  flame.  Even  though  it  is  fairly  constant  on  average,  /0.15 
exhibits  quite  strong  variability.  In  all  three  simulations  and  at  all  times  after  the  flame  has  reached  its  quasi-steady 
state,  /0.15  remains  above  unity  becoming  as  high  as  2.0  in  SI  and  1.5  in  S3.  Similarly  to  other  quantities  we  have 
considered,  /q.  15  shows  less  variability  with  increasing  resolution.  Time-averaged  values  of  /0.15  are  given  in  Table  3. 
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In  particular,  Z0.15  =  1.30  in  S3,  and  this  value  can  be  viewed  as  converged  to  within  a  few  percent  with  the  overall 
rate  of  convergence  being  somewhat  faster  than  linear. 

In  order  to  assess  the  dependence  of  /  on  the  choice  of  the  value  of  Y ,  consider  the  time-averaged  distributions 
of  I(Y)  shown  in  Fig.  4b.  In  particular,  we  show  I(Y)  =  (S  7 / S 7) / (At(Y) / L2).2  In  all  simulations,  7  is  substantially 
larger  than  unity  throughout  the  reaction  zone.  Also,  7  >  1.0  for  all  values  of  7,  with  the  only  exception  being  the 
outermost  region  of  the  preheat  zone  in  S3,  where  7  <  1.0  for  Y  >  0.93.  In  the  higher  resolution  simulations  S2  and 
S3,  7  decreases  monotonically  with  increasing  values  of  Y .  Similarly  to  A,  I  varies  the  least  inside  the  reaction  zone, 
i.e.,  for  Y  ~  0.15-0.5.  For  example,  the  difference  between  7o.is  and7o.s,  which  corresponds  to  the  boundary  between 
the  reaction  and  preheat  zones,  is  «2%  in  SI,  «5%  in  S2,  and  ^8%  in  S3. 

Finally,  the  ±cr  standard  deviation  of  individual  values  of  i(Y)  in  S3  is  shown  in  Fig.  4b  as  a  shaded  gray  area. 
The  <x  is  the  smallest  in  the  reaction  zone,  where  it  is  <  15%  of  7.  It  becomes  >  30%  in  the  preheat  zone  as  a  result  of 
the  stronger  variability  of  A  at  large  values  of  Y  (cf.  Fig.  3).  Note  that  throughout  the  reaction  zone,  i.e.,  for  Y  <  0.5, 
the  7  +  cr  interval  is  above  unity. 

3.3.  What  is  the  flame  surface  area? 

Distributions  of  the  time-averaged  stretch  factor  7,  shown  in  Fig.  4b,  fail  to  answer  the  question  concerning  the 
validity  of  Damkohler’s  concept  in  the  regime  considered  here.  Large  values  of  7,  found  in  the  reaction  zone,  indicate 
a  significant  departure  from  Damkohler’s  model.  At  the  same  time,  values  of  7  close  to  unity,  present  in  the  outer 
regions  of  the  preheat  zone,  suggest  that  the  increase  in  S  7  can  be  completely  accounted  for  by  the  increase  in  A. 

This  demonstrates  the  importance  of  the  question  previously  raised  in  §  1:  What  is  the  surface  area  of  the  flame 
in  such  high-speed  turbulent  flows?  In  particular,  since  the  response  of  different  parts  of  the  flamelet  to  the  action  of 
intense  turbulence  is  pronouncedly  distinct,  what  value  of  Y  most  fully  and  accurately  represents  the  overall  behavior 
of  the  turbulent  flame? 

These  questions  can  be  addressed  by  considering  the  correlation  of  A(Y )  with  two  global  properties  of  the  turbulent 
flame,  its  speed  and  width.  Fig.  5  shows  correlations  between  St/S 7  and  A(Y)/L2  calculated  for  three  representative 
values  of  Y:  Y  =  0.15,  corresponding  to  the  region  of  peak  reaction  rate,  Y  —  0.5,  corresponding  to  the  boundary 
between  the  reaction  and  the  preheat  zones,  and  Y  =  0.95,  corresponding  to  the  coldest  part  of  the  preheat  zone. 
The  left  panels  of  Fig.  5  show  the  time-evolution  of  all  quantities,  while  the  right  panels  show  the  corresponding 
correlation  scatter  plots. 

Since  both  the  local  flame  speed  and  the  induction  times  of  the  unburned  fuel  have  finite  values,  there  must  be  a 
delay  in  the  response  of  5  7  to  the  changes  in  the  flame  configuration,  represented  by  its  surface  area.  Therefore,  the 
degree  of  correlation  between  S  j  and  A(Y)  must  be  a  function  of  a  time  lag,  At,  between  these  two  quantities.  For  each 
value  of  Y  in  Fig.  5,  we  determined  the  time  lag,  At,  ,  that  produced  the  best  correlation.  This  was  done  by  directly 
calculating  the  cross-correlation  between  Sr(t)/S  7  and  A(Y,t)/L 2  and  finding  its  maximum.  The  value  obtained  for 
A tc  was  verified  by  determining  the  least  squares  fit  for  the  distribution  of  S  i  (t)/S  7  vs.  A(Y,t-  At')/ Lr  and  by  finding 
At'  that  maximized  the  slope  of  the  fit  and  minimized  its  residuals.  The  values  of  A(Y)/L 2  in  Fig.  5  were  then  shifted 
in  time  with  respect  to  St/S  l  by  the  corresponding  time  lag  A tc  given  in  the  lower  right  corner  of  the  panels  (a),  (c), 
and  (d). 

The  first  key  conclusion  emerging  from  Fig.  5  is  that  S  7  and  A  become  progressively  less  correlated  with  increas¬ 
ing  values  of  Y .  While  the  correlation  is  pronounced  throughout  the  reaction  zone  and  is  exceptionally  strong  near  the 
peak  reaction  rate,  S  7  and  A  are  only  very  weakly  correlated  in  the  colder  parts  of  the  preheat  zone.  Note  that  panels 
(b),  (d),  and  (f)  have  the  same  scale,  which  shows  the  relative  increase  in  the  scatter  of  the  distribution  and  its  overall 
shift  to  higher  values  of  A  in  the  preheat  zone. 

Note  in  Fig.  5b  that  at  the  peak  reaction  rate,  the  scatter  of  the  distribution  is  the  largest  near  the  average  values  of 
S  t  and  A0.15.  At  extreme  values,  however,  the  scatter  becomes  very  small  with  data  points  located  very  close  to  the 
fit  line.  This  effect  is  much  less  pronounced  for  Y  =  0.5,  and  for  Y  =  0.95  the  scatter  appears  to  increase  substantially 
toward  extreme  values  of  ,4 0.95. 

The  time  delays  that  give  the  best  correlation  between  the  two  quantities  are  positive,  and  thus  they  show  that  S  7 
does  indeed  lag  behind  A.  Moreover,  the  magnitude  of  A tc  increases  with  increasing  Y .  Qualitatively,  this  agrees  with 
the  fact  that  5  7  at  a  given  moment  is  primarily  determined  by  the  region  of  the  highest  reaction  rate,  which  corresponds 
to  the  lowest  values  of  Y .  Consequently,  the  time  delay  between  5  7  and  the  area  of  the  isosurface  representing  the 
peak  reaction  rate  is  very  small,  namely  0.04rej.  Since  turbulence  reorganizes  the  flame  on  a  timescale  ~  rec 7,  such  a 
small  lag  is  too  short  for  the  flame  structure  to  change  in  any  significant  way.  Therefore,  the  correlation  between  A0.15 
and  S  7  is  very  strong. 


-Note  that  using  instead  the  expression  I(Y)  =  (.S’ 7 IS  i)j(Aj( Y)j lr )  gives  virtually  the  same  result. 
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At  the  same  time,  since  the  reaction  rate  is  fairly  low  at  Y  —  0.5,  the  contribution  of  this  region  to  the  overall 
turbulent  flame  speed  is  small.  Burning  still  needs  to  accelerate  substantially  in  this  area  in  order  for  it  to  become  the 
new  region  of  peak  reaction  rate,  which  will  then  predominantly  determine  the  magnitude  of  5 7.  The  time  required 
for  this  to  occur  results  in  a  longer  time  delay,  namely  0.1  red.  During  this  time,  however,  turbulence  is  able  to  change 
the  flame  structure  more  strongly  than  in  the  case  of  Y  =  0.15.  As  a  result,  the  degree  of  correlation  between  A0.5  and 
5  t  decreases. 

In  the  coldest  regions  of  the  preheat  zone,  there  is  no  correlation  between  5  7  and  A0.95  at  small  values  of  A tc.  We 
found  evidence  of  weak  correlation  for  a  much  larger  time  delay  A tc  =  0.58rej.  This  value,  however,  is  comparable 
to  the  time  during  which  the  turbulence  completely  reorganizes  the  structure  of  the  turbulent  flame.  Therefore,  by 
the  time  burning  reaches  the  reactants  in  the  outer  regions  of  the  preheat  zone,  the  connection  between  their  original 
distribution  and  the  resulting  turbulent  flame  speed  is  significantly  disrupted. 

We  also  found  weak  correlation  between  A0.95  and  A0.15  with  the  time  lag  for  A0.95  similar  to  that  in  the  case  of 
S  7 ,  namely  Atc  =  0.56t„/,  and  the  slope  of  the  least  squares  fit  of  0.4907.  This  weak  correlation  between  /1 0.1 5  and 
A0.95  is  in  agreement  with  a  similarly  weak  correlation  between  5  7  and  A0.95.  In  particular,  it  reflects  the  fact  that 
after  the  time  ~  0.56  -  0.5Sred,  burning  reaches  the  outer  part  of  the  preheat  zone.  As  a  result,  the  former  Y  =  0.95 
region  becomes  the  new  Y  =  0.15  region  and  thus  the  new  site  of  the  peak  reaction  rate  which  now  determines  St- 
Substantial  change  in  the  flame-brush  structure  during  this  time,  however,  leads  only  to  weak  correlation  between 
A0.15  and  A0.95. 

The  computed  time  delays  can  be  compared  with  the  characteristic  propagation  time  of  the  laminar  flame, 


Tg 


Sl_ 

Sl 


30  L 

Yu 


=  3.75red, 


and  with  the  adiabatic  induction  time  for  the  reaction  model  used  in  this  work. 


'ind 


(13) 


(14) 


where  T,„d  is  based  on  the  Frank-Kamenetskii  approximation  [23].  In  eq.  (13),  we  used  the  values  of  L  and  U  given 
in  Table  2. 

Eq.  (14)  is  derived  under  the  assumption  that  heating  of  a  fluid  element  by  chemical  energy  release  greatly  ex¬ 
ceeds  its  cooling  due  to  thermal  conduction.  Therefore,  eq.  (14)  applies  only  in  the  reaction  zone.  There  Tmd  is  the 
controlling  timescale  and  it  is  1  -  2  orders  of  magnitude  smaller  than  Tg.  For  a  planar  laminar  flame. 


Tm  ~  O.OItj  »  0.05r«/,  for  Y  =  0.15, 
Tind  *  0.06r^  »  0.21r«/,  for  Y  =  0.5. 


(15) 


This  shows  that  for  Y  —  0.15,  Tlnd  is  practically  equal  to  the  corresponding  time  lag,  while  for  Y  =  0.5,  Tj„d  is  within 
a  factor  of  »  2  of  A tc  (Fig.  5a,c).  Using  in  eq.  (14)  the  specific  heat  at  constant  volume,  Cv  —  R/M(y  -  1),  instead  of 
Cp,  i.e.,  assuming  that  burning  occurs  at  constant  volume  rather  than  at  constant  pressure,  would  decrease  the  values 
of  Tind  by  ~  15%  to  0.04rfj  and  0.\?>Ted,  respectively. 

In  the  preheat  zone,  rmd  rapidly  becomes  much  larger  than  Tg,  since  heating  and  ignition  of  the  reactants  now 
are  controlled  by  heat  transfer  and  the  propagation  of  the  flame  as  a  whole.  For  Y  =  0.95,  however,  the  time  delay 
A tc  =  0.5 8 *  0. 1 6rn.  Furthermore,  it  is  an  order  of  magnitude  smaller  than  the  time  r  »  1 ,676l/S  l  =  1 .67  Tg, 
which  is  necessary  for  the  flame  to  propagate  over  the  characteristic  distance  separating  Y  =  0.95  and  Y  =  0. 1 5  in  the 
flamelet  structure  (cf.  Fig.  7  in  [9]).  This  shows  that  while  in  the  reaction  zone  the  time  delays  are  in  good  agreement 
with  the  corresponding  Tmd,  in  the  preheat  zone  A tc  is  substantially  smaller  than  the  characteristic  propagation  times 
of  the  laminar  flame. 

In  addition  to  the  strong  correlation  between  the  values  of  5  7  and  A0.15,  the  distribution  shown  in  Fig.  5b  has 
another  important  property  that  is  not  present  at  higher  values  of  Y.  As  the  flame  becomes  less  convolved,  it  will 
behave  progressively  more  like  a  planar  laminar  flame.  This  means  that  Ay/L 2  for  the  isosurface  most  representative 
of  the  flame  surface  and  S  7  /S  7  must  approach  unity  simultaneously.  Consider  now  the  linear  least  squares  fit  for  the 
case  Y  =  0.15,  shown  with  a  solid  line  in  Fig.  5b.  It  has  the  form 


fr =  L44Tf 


0.41. 


(16) 


For  A0.15/X2  =  1,  this  expression  gives  S  7  /S  7  =  1.03,  which  is  the  deviation  of  only  3%  from  unity.  In  fact,  the  line 
which  passes  exactly  through  the  point  S  7  /S  /  =  A0.15/L2  =  1  and  also,  as  required  for  the  least  squares  fit,  through 


the  center  of  mass  of  the  distribution,  i.e.,  the  point  St /Sl  =St/Sl  andAo.is/L2  =  Ao.i5/L2,is 
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1  (A0.15/L2)  -  1 

=  1.46^-0.46. 


(17) 


Here  St/S y  and  A q  \s/L2  are  the  time-averaged  values  given  in  Table  3.  Expressions  for  the  coefficients  in  eq.  (17) 
approximate  the  first  and  second  coefficients  in  eq.  (16)  with  the  accuracy  of  w  1%  and  ~  12%,  respectively. 

This  demonstrates  that  the  distribution  of  S t/S l  as  a  function  of  A0.15/L2  (Fig.  5b)  closely  follows  the  trend  that 
has  the  correct  limiting  behavior  as  A0.15/L2  — >  1.  Linear  least  squares  fits  for  the  two  other  values  of  Y,  shown  in 
Fig.  5d,f,  do  not  recover  the  value  of  S  y  /.S'  /  =  1  as  A(Y)/L?  — >  1 .  It  is  particularly  important  to  consider  such  limiting 
behavior  for  Y  =  0.5,  where  a  reasonably  good  correlation  between  S  7  and  /to. .5  could  suggest  that  this  value  of  Y  can 
also  be  viewed  as  representative  of  the  flame  properties. 

Next  consider  the  correlation  of  A(Y )  with  the  second  key  global  characteristic  of  the  turbulent  flame,  namely  its 
width,  6t .  Fig.  6  shows  the  correlation  between  A(Y)/Lr  and  by  j L  for  the  same  three  values  of  Y  as  in  Fig.  5.  Since 
both  A  and  5t  represent  the  same  instantaneous  configuration  of  the  flame,  no  time  delay  between  them  would  be 
expected.3  Consequently,  no  time  shift  was  applied  to  either  quantity  in  Fig.  6. 

The  distribution  of  dy/L  as  a  function  of  A(Y)/L2  shows  the  same  trend  observed  for  S  y  /S  / .  In  particular,  the 
correlation  between  the  two  quantities  is  also  the  strongest  for  Y  =  0.15,  and  it  decreases  with  increasing  Y .  The 
values  of  by  / L  and  A0.95  appear  to  be  almost  completely  uncorrelated.  Note,  however,  that  overall  the  correlation  in 
the  reaction  zone,  and,  in  particular,  in  the  region  of  peak  reaction  rate,  is  weaker  for  by  than  S  y  with  the  distributions 
of  by  having  a  much  larger  scatter. 

The  decrease  in  correlation  between  by/L  and  A{Y)/Lr  with  increasing  Y  reflects  the  fact  that  the  isosurfaces  of 
higher  Y  tend  to  be  wrinkled  on  progressively  smaller  scales,  as  discussed  in  §  3. 1  (also  see  [9]).  As  a  result,  while  the 
areas  of  these  isosurfaces  increase,  they  contribute  less  to  by  which  is  determined  predominantly  by  the  flame  folding 
on  the  largest  scales.  Therefore,  by  is  best  correlated  with  A{Y)  at  the  lowest  values  of  Y ,  where  the  folding  on  such 
large  scales  is  most  pronounced. 

The  strong  correlation  of  A0.15  with  5 y  and  by,  which  represent  both  the  global  energetics  and  the  global  structure 
of  the  turbulent  flame,  as  well  as  the  correct  limiting  behavior  of  the  distribution  of  St/S  l  as  A0.15/L2  — >  1  present  a 
compelling  case.  These  results  show  that  Y  =  0.15  most  accurately  characterizes  the  overall  evolution  of  the  flame. 
Consequently,  the  flame  surface  area,  Ay,  is  best  represented  by  the  isosurfaces  of  Y  close  to  the  peak  reaction  rate. 

3.4.  Stretch  factor  and  the  balance  between  S  y  and  At 

We  can  now  revisit  the  question  of  the  actual  values  of  the  stretch  factor,  /,  in  the  regime  considered  here.4  Our 
results  demonstrate  that  /0.15  must  be  viewed  as  most  characteristic  of  the  system  evolution.  The  time  history  of  /0.15 
(Fig.  4a)  and  its  time-averaged  values  (Table  3)  show  pronounced  deviation  from  Damkohler’s  concept  on  the  order 
of  30%  as  the  turbulent  flame  speed  exhibits  a  strong  exaggerated  response  to  the  increase  in  the  flame  surface  area. 

This  exaggerated  response  is  illustrated  in  Fig.  5  as  the  shaded  gray  area.  In  particular,  the  shaded  region  in  Fig.  5b 
is  a  measure  of  the  relative  increase  of  the  actual  values  of  St/S  t  with  respect  to  the  values  given  by  Ao.  15  /Lr ,  which 
are  shown  as  a  dash-dot  line.  At  larger  values  of  A0.15/L2,  the  overall  distribution  deviates  progressively  more  from 
the  St/S l  =  A0.15/L2  line.  This  demonstrates  the  following  crucial  fact:  the  magnitude  of  the  exaggerated  response 
of  St  grows  with  increasing  flame  surface  area.  We  examine  this  result  in  more  detail  in  §  4.3. 

Finally,  the  time-averaged  value  of  the  stretch  factor  7  as  1  in  the  coldest  region  of  the  preheat  zone  (Fig.  4b)  was, 
in  fact,  simply  the  result  of  time  averaging.  Fig.  5e,f  shows  that  /0.95  has  a  much  larger  range  of  variability  than  /0.15, 
namely  0.5  -  1.7  (also  consider  the  ±cr  standard  deviation  of  HY )  shown  in  Fig.  4b).  At  the  same  time,  in  a  system 
in  which  St  is  determined  only  by  At,  I  by  definition  must  be  unity  and  constant.  This  further  shows  that  7  ~  I  at 
large  values  of  Y  in  simulation  S3  cannot  serve  as  evidence  of  the  fact  that  the  observed  S  y  is  determined  only  by  the 
increase  in  the  flame  surface  area. 


■'As  a  check,  we  verified  that,  indeed,  the  cross-correlation  had  its  maximum  at  the  zero  time  lag. 

4In  light  of  the  fact  that  S  y  and  A(Y)  are  best  correlated  when  A  is  shifted  in  time  with  respect  to  Sr  -  the  question  arises  whether  such  time  shift 
should  also  be  applied  to  A  when  calculating  I.  i.e.,  whether  the  expression  t(Y.  t)  =  (SrU)lS  i)f(A[  Y.  t  -  Atc)/L2)  should  be  used.  In  principle,  an 
argument  could  be  made  that,  due  to  the  inherent  delay  in  the  response  of  S  y  to  the  changes  in  the  flame  structure,  the  current  value  of  S  y  is  the 
result  of  A  that  existed  earlier  in  time  and,  thus,  /  must  be  calculated  using  this  earlier  value.  At  the  same  time,  as  was  shown,  the  time  lag  A tc  is 
extremely  small  for  Y  =  0.15,  which  as  we  have  concluded  should  be  used  to  estimate  I.  During  this  time,  the  change  in  the  flame  surface  area  is 
negligible  and  applying  the  time  shift  in  the  expression  for  /  does  not  cause  any  appreciable  change  in  the  result. 
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4.  Flame  surface  density  and  its  relation  to  the  turbulent  flame  speed 


4.1 .  Surface  density  of  the  fuel  mass  fraction  isosurfaces 


The  analysis  of  the  flame  surface  area  does  not  address  the  question  of  how  the  flamelets  are  organized  inside  the 
flame  brush  and,  in  particular,  how  tightly  they  are  packed.  It  can  be  seen  in  Fig.  6b  that  the  relative  variation  of  A0.15 
is  «  50%  larger  than  that  of  by  with  the  ratio  of  the  maximum  and  minimum  values  of  ,4 0.1 5  being  =  3  as  opposed  to 
«  2  for  6t-  This  suggests  that  the  change  in  the  flame  surface  area  cannot  be  accounted  for  only  by  the  variations  in 
the  turbulent  flame  thickness,  but  it  must  also  be  related  to  how  the  flame  surface  is  folded  inside  the  flame  brush. 

Consider  the  surface  density  of  the  fuel  mass  fraction  isosurfaces,  defined  as  the  isosurface  area  normalized  by  the 
effective  flame-brush  volume 


£00  = 


MY) 
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Detailed  discussion  of  this  definition  of  X(K )  and,  in  particular,  of  the  choice  of  the  uniform  normalization  by  the  full 
flame-brush  volume  is  given  in  Appendix  A. 

The  evolution  of  X0.01  and  X0.99,  where  ’Ey  =  X(K  =  Y'),  is  shown  in  Fig.  7a, b.  The  overall  trends  are  similar  to 
those  exhibited  by  /to. 01  and  A0.99  (cf.  Fig.  3a, b).  In  particular,  X  is  lower  on  the  product  side  of  the  flame.  Moreover, 
similar  to  A,  X  appears  to  converge  at  lower  values  of  Y  and  shows  virtually  no  convergence  in  the  preheat  zone.  Most 
important,  there  is  substantial  variability  in  X,  especially  on  the  fuel  side,  where  X  changes  by  up  to  a  factor  of  3. 
This  tighter  packing  of  the  isosurfaces,  which  periodically  occurs  in  the  course  of  system  evolution,  suggests  that  the 
increase  in  the  flame  surface  area  is  indeed  associated  with  the  increase  not  only  of  6t,  but  also  of  X. 

Fig.  7c  shows  time-averaged  distributions  of  X(T)  in  simulations  SI  -  S3.  Comparison  with  Fig.  3c  demonstrates 
that  both  A(Y)/L 2  and  X(T)  show  remarkably  similar  behavior.  The  time-averaged  values  of  X  are  also  not  constant, 
but,  instead,  they  progressively  increase  through  the  flamelet  from  its  product  side  to  the  fuel  side.  The  resulting 
inverted-S  shape  of  the  profiles  is  very  similar  to  that  of  A(Y)  with  both  A  and  X  at  the  two  extreme  values  of  Y 
(Y  —  0.01  and  0.99)  differing  in  S3  by  almost  a  factor  of  two.  Similar  to  A,  X  varies  the  least  in  the  interior  of 
the  reaction  zone,  i.e.,  for  Y  w  0.15  -  0.6,  while  outside  this  region  it  shows  strong  dependence  on  Y .  Just  outside 
the  reaction  zone,  individual  profiles  begin  to  diverge  from  each  other  with  the  variation  between  them  becoming 
progressively  larger  for  higher  values  of  Y . 

In  the  reaction  zone,  values  of  X  for  all  three  calculations  are  much  closer  than  those  of  A.  In  fact,  Table  3  shows 
that  X(T)  exhibits  4,;' -order  convergence  compared  with  3rd-order  convergence  for  A.  As  a  result,  values  of  X  are 
virtually  identical  for  Y  ~  0.01  -  0.4  in  S2  and  S3.  On  the  other  hand,  in  the  preheat  zone,  the  profiles  of  X  diverge 
more  than  in  the  case  of  A.  This  shows  that  the  effects  of  small-scale  turbulent  motions  on  the  reaction  and  the  preheat 
zones  are,  in  fact,  even  more  disparate  in  terms  of  their  ability  to  provide  tight  folding  of  the  isosurfaces  compared 
with  their  ability  to  increase  the  isosurface  area. 

Flame  surface  density  is  a  quantity  used  both  in  experimental  and  theoretical  combustion  research.  In  particular, 
X„MV  is  defined  as  the  surface  density  in  the  region  of  the  flame  brush  with  mean  reactedness  c  -  0.5,  where  c  = 
(T  -  Tq)/(Tp  -  To).  It  can  then  be  shown  that  [3] 


X 
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where  /\'  5  is  the  flamelet  surface  area  based  on  c  -  0.5,  and  A  *  represents  the  average  area  of  the  flame  brush.  In 
a  system  described  by  the  one-step  Arrhenius  kinetics,  c  =  0.5  corresponds  to  Y  =  0.5.  Therefore,  A*5  =  A0.5. 
At  the  same  time,  in  a  planar  brush  A*L  =  L 2,  as  was  discussed  in  §  3.1.  Thus,  Xmav  is  equivalent  to  X0.5  as  given 
by  eq.  (18).  It  is  reasonable  then  to  compare  values  of  X0.5  in  our  simulations  with  those  of  X„MV  obtained  in  other 
experimental  and  numerical  settings.  In  particular,  in  S3,  X0.5  =  0.73  mm  This  is  comparable  to  the  range  of  values 
of  Emax  «  0.12-  0.6  mm~ 1  obtained  in  a  number  of  experimental  studies  using  a  wide  variety  of  flame  configurations 
(see  [3]  for  a  summary  table).  The  fact  that  the  value  in  S3  is  somewhat  larger  is  not  surprising,  given  that  the  turbulent 
intensity  in  it  is  substantially  higher. 

Overall,  however,  the  traditional  definition  of  X  given  in  [3]  is  not  equivalent  to  the  definition  used  here  unless 
Y  =  0.5.  Therefore,  unlike  X„MA,  X0.5  is  not  necessarily  the  maximum  value  of  X  for  all  values  of  Y .  Fig.  7c  shows  that 
while  in  SI  X0.5  is  indeed  approximately  the  largest  value,  this  is  not  the  case  in  S2  and  S3.  At  the  same  time,  in  all 
three  calculations,  X0.5  is  representative  of  the  values  of  X  throughout  the  reaction  zone  with  the  difference  between 
X0.5  and  X0.15  being  =6-9%  in  the  higher-resolution  cases. 


4.2.  Distributions  ofA(Y)  andE(Y)  and  the  effects  of  small-scale  turbulence 

How  do  these  results  concerning  the  distributions  of  A(Y )  and  X(K )  relate  to  previous  conclusions  [9]  regarding  the 
inability  of  the  turbulent  cascade  to  penetrate  the  flamelet  interior  and  the  role  of  small-scale  turbulent  motions?  As 
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was  discussed  in  §  2. 1  (also  see  [9]  for  further  details),  the  progressive  increase  in  resolution  from  SI  to  S3  extends  the 
turbulent  cascade  to  smaller  scales,  and  this  leads  to  a  substantial  increase  in  the  energy  of  turbulent  motions  on  scales 
A  <  6l  in  nonreactive  turbulence.  This  increase  in  resolution  was  found  to  cause  the  flame  surface  on  the  fuel  side  to 
be  wrinkled  on  progressively  finer  scales,  while  remaining  virtually  unchanged  on  the  product  side.  It  was  determined 
that  the  deviation  of  the  internal  flamelet  structure  from  that  of  the  planar  laminar  flame  increases  with  decreasing 
temperature.  These  two  results  showed  that  the  effects  of  small-scale  turbulent  motions  are  most  pronounced  in  the 
coldest  parts  of  the  preheat  zone.  These  effects  diminish  with  increasing  temperature  and  completely  disappear  once 
the  reaction  rate  becomes  significant  [9]. 

Results  presented  in  this  paper  are  consistent  with  this  picture.  The  steady  growth  of  A/L 2  and  X  on  the  fuel  side 
of  the  preheat  zone  with  increasing  resolution  shows  that  more  intense  small-scale  motions  indeed  fold  isosurfaces 
on  progressively  finer  scales.  It  must  be  emphasized  that  the  change  in  A  alone  does  not  necessarily  mean  that  it  is 
caused  by  finer  wrinkling  on  smaller  scales.  Only  when  viewed  in  conjunction  with  the  surface  density,  which  shows 
a  very  similar  dependence  on  Y,  does  such  increase  in  A  serve  as  a  strong  indication  of  the  small-scale  wrinkling. 

With  increasing  temperature,  not  only  do  A/L 2  and  X  decrease  rapidly,  but  also,  most  importantly,  individual 
profiles  for  S2  and  S3  begin  to  approach  each  other.  This  demonstrates  that  the  main  difference  between  these  two 
simulations,  namely  the  presence  of  more  energetic  small-scale  turbulence  in  S3,  is  being  eliminated  and,  thereby, 
motions  on  scales  A  <  0/  are  gradually  suppressed.  Since  it  is  these  small  scales  that  enhance  the  diffusive  transport, 
which  in  turn  broadens  the  preheat  zone,  their  suppression  causes  the  internal  flame  structure  to  approach  that  of  the 
planar  laminar  flame  [9]. 

As  the  reaction  rate  becomes  substantial  at  Y  <  0.5,  both  A/L2  and  X  become  similar  in  S2  and  S3.  At  this  point  the 
only  energetic  scales  are  A  >  6l,  which  are  originally  the  same  in  both  calculations  (cf.  Fig.  1  in  [9]).  Consequently, 
they  generate  similar  isosurface  areas  and  densities.  Since  these  scales  cannot  support  small-scale  diffusive  transport, 
any  broadening  of  the  internal  flamelet  structure  effectively  disappears  at  Y  <  0.5  [9].  Further  decrease  in  A/L2  and 
X  is  the  consequence  of  the  continued  suppression  of  the  progressively  larger  scales  A  >  6l- 

There  is  one  important  distinction  between  the  distributions  of  A/L2  and  X  and  the  time-averaged  profiles  of  Y 
and  T,  which  represent  the  internal  flamelet  structure  [9],  The  Y  and  7  profiles  showed  very  little  variation  between 
simulations  SI -S3  not  only  in  the  reaction  zone,  but  also  in  the  preheat  zone.  This  is  in  contrast  with  the  distributions 
A/L 2  and  X,  which  differ  substantially  in  the  preheat  zone  even  between  the  S2  and  S3. 

This  suggests  that  small  scales,  A  <  6l,  contribute  differently  to  the  turbulent  diffusive  transport  and  to  the  wrin¬ 
kling  of  the  isosurfaces.  The  former  is  primarily  governed  by  the  largest  of  these  small  scales,  i.e.,  scales  not  much 
smaller  than  b/,  since  these  scales  have  the  highest  velocities  associated  with  them.  The  energy  on  these  scales  is  the 
same  in  all  three  simulations  and,  thus,  they  produce  the  same  structure  of  the  broadened  preheat  zone.  At  the  same 
time,  all  scales  smaller  than  6l  contribute  to  the  isosurface  wrinkling  and  the  different  energy  content  of  these  scales 
causes  the  distributions  of  A/L 2  and  X  to  differ  substantially  in  the  preheat  zone. 

Qualitatively,  both  the  internal  flamelet  structure  described  in  [9]  and  the  distributions  of  A  and  X  presented  here 
create  a  consistent  picture  of  the  profound  transformation  that  turbulence  undergoes  as  it  passes  through  the  flame. 
Changes  in  the  energy  budget  between  different  scales,  which  are  most  certainly  accompanied  by  the  development  of 
both  anisotropy  and  inhomogeneity  of  the  velocity  field,  are  complex.  While  the  evidence  presented  here  and  in  [9] 
provide  hints  about  the  nature  of  this  transformation,  the  details  remain  unknown.  For  instance,  are  the  energy  release 
and  the  resulting  fluid  expansion  the  only  effects  responsible  for  altering  the  turbulent  field  and  redistributing  energy 
between  different  scales?  What  are  the  relative  contributions  of  various  small  scales  to  the  increase  in  A  and  X?  How 
does  the  shift  in  balance  between  small  and  large  scales  with  increasing  temperature  change  at  different  turbulent 
intensities?  All  of  these  questions  need  to  be  addressed  in  future  studies. 

4.3.  Correlation  between  S t  and  X(T) 

Peaks  and  troughs  in  values  of  Xo.oi  in  Fig.  7a  can  generally  be  associated  with  those  in  values  of  5  j-  in  Fig.  lb.  This 
suggests  that  the  density  of  the  flame  surface,  packed  inside  the  flame  brush,  is  a  dynamically  important  characteristic 
affecting  the  energetics  of  the  system  and,  thus,  its  evolution.  In  order  to  consider  this  question  more  rigorously,  we 
show  in  Fig.  8  the  correlation  between  St/S l  and  X(T)  for  the  same  three  values  of  Y  as  in  Figs.  5-6. 

As  with  A,  one  would  expect  St  to  respond  with  a  certain  delay  to  the  changes  in  X.  Even  though  both  A  and  X 
represent  the  same  instantaneous  flame  structure,  we  chose  not  to  make  an  a  priori  assumption  that  the  magnitude  of 
that  delay  is  the  same  as  it  was  in  the  case  of  A.  We  determined  the  time  lag,  A  tc,  which  produced  the  best  correlation 
between  5  T(t)/S l  and  X(K,  t),  for  each  individual  value  of  Y  following  the  same  procedure  as  was  used  in  §  3.3.  The 
values  of  X(T)  in  Fig.  8  were  shifted  in  time  by  At,  given  in  the  lower  right  corner  of  the  panels  (a),  (c),  and  (d). 

The  principal  result  is  the  strong  correlation  of  the  values  of  X0.15,  corresponding  to  the  region  of  the  peak  reaction 
rate  in  the  flamelet,  with  S  t  ■  The  scatter  of  the  distribution  of  S  ■/■  /S  /  as  a  function  of  X0.15,  however,  is  larger  than  for 
/1 015.  At  Y  —  0.5,  the  correlation  of  A  7  with  X0.5  is  weaker,  as  for  A0.5,  although  it  is  still  quite  pronounced.  Finally, 
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similar  to  A0.95 ,  £0.95  shows  at  most  only  very  weak  correlation  with  Sy  at  A tc  —  0.58 Tej  (Fig.  8f).  In  addition,  we 
found  some  evidence  of  weak  anticorrelation  of  E0.95  with  .S'  7  at  At,  -  0,  with  the  slope  of  the  least  squares  fit  of 
-3.007.  This  is  in  contrast  to  A0.95,  which  was  completely  uncorrelated  with  St  at  values  of  At,  close  to  zero.  It 
is  not  clear,  however,  what  the  physical  significance  of  this  anticorrelation  is,  since  it  is  unlikely  for  5  t  to  have  an 
instantaneous  response  to  the  distribution  of  cold  unreacted  fuel.5 

The  time  delays  for  2  in  the  reaction  zone  are  similar  to  those  obtained  in  the  case  of  A.  Moreover,  A tc  for  E0.5 
is  much  closer  to  the  magnitude  of  the  induction  time  given  in  eq.  (14),  than  it  was  for  A0.5,  being  within  25%  of  the 
corresponding  value  of  Tmci- 

These  results  are  further  evidence  that  Y  close  to  the  peak  reaction  rate  is  most  representative  of  the  dynamics  of 
the  turbulent  flame.  Therefore,  as  A0.15  was  concluded  in  §  3.3  to  represent  the  flame  surface  area,  X0.15  can  be  viewed 
to  represent  the  flame  surface  density,  Ey .  Consequently,  the  fact  that  larger  values  of  S  y  are  correlated  with  larger 
values  of  X0.15  shows  that  the  increase  of.S'  y  is  closely  associated  with  more  tightly  packed  flame  configurations. 

So  far  we  have  considered  the  correlation  between  St,  A0.15,  and  E0.15  separately.  In  order  to  demonstrate  the 
close  connection  between  these  three  quantities,  we  colored  each  data  point  in  Fig.  3b  according  to  the  corresponding 
magnitude  of  2o.  15 .  It  can  be  seen  that,  indeed,  larger  values  of  S  t  and  Ao.  15  are  associated  with  larger  values  of  2o.  15 . 

In  particular,  an  increase  in  the  value  of  /0.15  from  ~  1.2  to  ~  1.36  is  associated  with  the  increase  in  X0.15  by  more  than 
a  factor  of  two  from  x;0.4  mm-1  to  r;0.85  mm-1  . 

This  finally  leads  us  to  the  following  two  key  conclusions: 

1 .  Larger  values  of  At  are  associated  with  larger  values  of  Ey .  This  shows  that  increase  in  the  flame  surface  area 
in  the  course  of  the  flame  evolution  is  primarily  due  to  the  much  tighter  folding  and  packing  of  the  flame,  rather 
than  only  due  to  the  increase  in  the  overall  extent  of  the  flame  brush. 

2.  Increase  in  At  and  Ey  leads  not  only  to  larger  values  of  S  y  but,  most  importantly,  also  to  a  larger  deviation  of 
St/S l  from  Ay/L2,  i.e.,  from  the  value  predicted  by  Damkohler’s  concept.  In  other  words,  in  the  presence  of 
the  high-speed  turbulence,  the  increase  in  the  flame  surface  area  and  the  associated  tighter  packing  of  the  flame 
result  in  the  progressively  more  exaggerated  response  of  St  to  such  increase  in  Ay-  This,  therefore,  causes  a 
substantially  accelerated  burning. 

5.  Mechanism  of  the  turbulent  flame  speed  increase 

5.1 .  Potential  causes  of  the  increase  of  the  local  burning  speed 

The  substantial  deviation  of  1  from  unity  observed  in  the  simulations  brings  up  the  following  question:  What 
mechanism  is  responsible  for  raising  5  j  beyond  what  can  be  attributed  to  the  increase  in  Ay? 

According  to  Damkohler’s  concept,  I  >  1  results  from  the  change  in  the  local  flame  speed,  S 7,  which  is  the 
consequence  of  the  enhancement  of  diffusive  processes  by  turbulent  transport.  Therefore,  this  change  cannot  be 
episodic  and  local,  but  rather  it  must  be  a  reflection  of  the  statistically  dominant  state.  It  was  demonstrated  in  [9], 
however,  that  there  is  no  evidence  of  such  enhancement,  as  the  average  internal  flamelet  structure  in  the  reaction  zone 
is  virtually  identical  to  that  of  the  planar  laminar  flame,  and  the  resulting  thermal  width  of  the  flamelets  is  practically 
equal  to  by.  Since  for  the  first-order  Arrhenius  kinetics  there  is  a  unique  correspondence  between  the  flame  structure 
and  speed,  this  showed  that,  on  average,  the  flame  propagates  locally  with  the  speed  close  to  Sy.  Furthermore,  there 
is  no  significant  fuel  preconditioning  by  the  turbulent  flow,  which  could  increase  S  7.  For  the  full  duration  of  the 
simulations,  the  pressure  in  the  domain  remains  smooth  and  close  to  its  initial  value  Pq.  The  fuel  temperature  rises  by 
the  end  of  the  simulation  to  ^330K,  which  is  not  sufficient  to  increase  5  7  appreciably. 

These  considerations  show  that  the  observed  anomalous  increase  in  5  j  cannot  be  accommodated  within  Damkohler’s 
paradigm.  Consequently,  a  different  mechanism  must  augment  the  effect  of  the  increased  flame  surface  area. 

Instead  of  being  the  result  of  a  uniform  global  increase  of  5  7  compared  to  5  y  throughout  the  flame  brush,  enhanced 
values  of  5  y  could  be  caused  by  a  significant  local  increase  in  5  7.  Two  immediate  causes  for  such  local  enhancement 
can  be  envisaged:  (1)  local  inhomogeneities  in  the  thermodynamic  state  of  the  fuel,  and  (2)  intermittency  of  the 
turbulent  flow. 

We  did  not  find,  however,  either  of  them  to  be  a  significant  factor  in  causing  the  observed  increase  in  Ay.  In 
particular,  in  a  subsonic  turbulent  flow,  such  as  the  one  considered  here,  with  the  turbulent  Mach  number  Ma  substan¬ 
tially  less  than  unity  (see  Table  1),  local  variations  of  density  or  temperature  are  much  too  small  to  affect  S  7  in  any 
appreciable  way. 


5We  note  that  we  were  unable  to  find  any  statistically  significant  correlation  or  anticorrelation  between  So .15  and  S095  at  any  time  lag.  This  is 
in  contrast  with  the  weak  correlation  between  A0.15  and  A0.95  discussed  in  §  3.3. 
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The  potential  role  of  turbulence  intermittency  on  the  flame  has  been  studied  by  Pan  et  al.  [24]  in  the  context 
of  the  Rayleigh-Taylor-driven  thermonuclear  flames  in  the  interior  of  a  white  dwarf  during  the  Type  la  supernova 
explosion.  Their  results  suggested  the  possibility  that  intermittency  in  the  velocity  field  causes  local  disruption  and 
broadening  of  the  flame,  which  potentially  would  also  increase  the  local  burning  speed.  It  is  not  immediately  clear, 
however,  to  what  degree  their  conclusions  would  apply  to  the  system  discussed  here.  The  conditions  in  the  white 
dwarf  interior  considered  in  [24]  are  substantially  different  from  the  ones  in  SI  -  S3.  Moreover,  their  analysis  did  not 
consider  any  feedback  of  the  flame  on  the  turbulent  field  and,  in  particular,  it  did  not  include  fluid  expansion  due  to 
heat  release  which  invariably  affects  intermittency  of  the  flow.  While  we  do  observe  intermittency  in  the  flow  field  in 
our  simulations,  we  find  that  regions  of  large  velocity  enhancement  are  statistically  too  rare,  their  spatial  extent  is  too 
small,  and  they  are  much  too  short  lived  to  increase  S  /  and,  thus,  to  make  any  significant  contribution  to  the  overall 
increase  in  Sr-  At  the  same  time,  the  role  of  intermittency  in  the  dynamics  of  turbulent  flames  does  require  further 
investigation.  In  particular,  it  would  be  extremely  important  to  extend  the  analysis  of  Pan  et  al.  [24]  by  including  the 
feedback  of  the  flame  on  turbulence  under  conditions  characteristic  of  chemical  rather  than  thermonuclear  flames. 

5.2.  Local  increase  ofS  /  in  cusps  formed  by  flame  collisions 

According  to  the  theory  of  flame  stretch  (see  [5]  for  a  review),  under  the  Le  -  1  conditions,  the  flame  is  not 
affected  either  by  the  flow-induced  strain  or  curvature,  in  accordance  with  Damkohler’s  concept.  In  particular,  a 
stationary  spherical  flame  supported  by  a  point  source  of  mass,  which  is  an  example  of  a  curved  unstrained  flame,  has 
the  internal  structure,  and  thus  local  burning  velocity,  identical  to  that  of  a  planar  laminar  flame  [5],  Such  analysis, 
however,  is  typically  performed  in  the  regime  when  curvature  radius  rc  »  5/ .  Our  simulations  of  an  idealized 
spherical  or  cylindrical  flame  collapsing  into  a  stationary  fuel  corroborate  this  result.  Furthermore,  they  show  that 
Si  =  Sl  not  only  when  rc  »  5/,,  but  up  until  the  moment  when  the  flame  curvature  becomes  «  I /c)/,,  i.e.,  until 
the  flame  effectively  collapses  onto  itself.  At  this  point,  S  /  increases  substantially.  Beyond  this  moment,  however, 
the  flame  itself  effectively  ceases  to  exist  and,  therefore,  the  local  flame  speed  looses  its  meaning.  Such  idealized 
situation  is  hardly  representative  of  the  conditions  arising  in  an  actual  turbulent  flame.  Nonetheless,  regions  of  large 
flame  curvature  ~  1/5/,  are  frequently  created  in  the  high-speed  turbulent  flow  considered  here,  and  we  now  show  how 
they  naturally  provide  a  mechanism  for  a  significant  local  enhancement  of  S  /. 

We  have  demonstrated  above  that  in  the  presence  of  intense  turbulence,  larger  values  of  I  are  well  correlated  with 
the  increase  in  X  and,  thus,  with  much  tighter  folding  and  packing  of  the  flamelets  inside  the  flame  brush.  The  inverse 
of  the  surface  density  is  a  measure  of  an  average  separation  between  surface  elements.  Given  that  Xo.15  =  0.67  mm-1 
in  simulation  S3,  the  average  distance  between  individual  flame  sheets  is  l/Xo.15  ~  1-49  mm  ~  4.7c)/,.  At  the  same 
time,  X0.15  can  be  as  high  as  0.97  mm'1  (Figs.  5b,  8b),  resulting  in  an  even  lower  average  separation  of  »  3.2c)/,.  This 
is  comparable  to  the  full  width  of  individual  flamelets  (cf.  Fig.  7  in  [9]).  Such  tightly  folded  flamelet  configurations 
invariably  result  in  frequent  collisions  of  individual  flame  sheets. 

Fig.  9  gives  an  example  of  such  collision  in  simulation  S3.  The  figure  shows  the  flame-brush  structure  at  three 
times:  1 1.86 Tej  (upper  panel),  as  well  as  0.1  Ted  —  3  ps  and  0.2 Ted  =  6  ps  later  (middle  and  lower  panels,  respectively). 
Initially,  a  region  with  highly  convolved  flame  develops  in  the  flame  brush  (region  A).  At  this  time,  S 1  =  S /,  locally 
throughout  region  A  and  self-propagation  of  the  flame  can  be  neglected  on  timescales  considered  in  Fig.  9.  As 
turbulent  motions  continue  to  bring  individual  flame  sheets  closer  to  each  other,  the  curvature  radius  of  the  flame 
becomes  close  to  5/ ,  and  the  preheat  zones  begin  to  overlap  substantially  over  an  extended  region  of  the  flame  surface. 
This  marks  the  formation  of  two  regions  of  high  flame  curvature  >  1  /5/,  (regions  B),  which  we  will  hereafter  refer  to 
as  “cusps.”  Regions  C  in  the  lower  panel  show  that  3  ps  later  the  flame  sheets  have  merged  and  formed  two  extended 
reaction  zones,  which  suggests  substantially  accelerated  burning  in  that  area  of  the  flame  brush.  Note  also  a  rapid 
decrease  of  the  area  of  the  Y  =  0.5  isosurface  represented  with  the  thin  black  line.  Comparison  of  regions  B  and  C 
shows  that  this  isosurface  propagated  over  the  distance  ~  0.5  mm  over  the  time  3  ps,  which  implies  the  propagation 
speed  (but  not  the  local  burning  speed)  «55 S  l.  During  the  time  shown  in  Fig.  9,  two  other  highly  elongated  regions 
of  flame  collision  have  formed  near  the  upper  face  of  the  domain.  Thus,  Fig.  9  illustrates  the  fact  that  flame  collisions 
are  ubiquitous  in  the  turbulent  flame,  and  they  do  produce  regions  of  large  flame  curvature.  Moreover,  their  evolution 
in  Fig.  9  suggests  substantial  increase  of  both  the  local  burning  and  propagation  speeds  in  the  cusps. 

5.3.  Structure  and  properties  of  cusps 

Cusps  were  first  considered  by  Zel’dovich  [25]  (also  see  [26]),  who  suggested  them  as  a  mechanism  of  flame 
stabilization  that  prevents  the  unbounded  exponential  growth  of  the  flame  surface  under  the  action  of  the  Landau- 
Darrieus  instability  [27].  Such  regions  would  have  a  high  propagation  velocity,  causing  rapid  annealing  of  the  flame 
surface,  and  thus  providing  an  efficient  mechanism  to  counterbalance  the  flame  instability  [25],  It  was  also  suggested 
in  [25]  that  S  /  near  the  tip  of  the  cusp  ceases  to  be  equal  to  S  /,,  and  the  cusp  region  develops  a  structure  similar  to 
the  tip  of  the  Bunsen  flame.  The  flame  in  the  analysis  given  in  [25,  26],  however,  was  considered  to  be  a  gasdynamic 
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discontinuity,  and  so  its  internal  structure  was  ignored.  As  a  result,  the  actual  increase  of  the  local  flame  speed  in  the 
cusp  could  not  be  determined. 

In  order  to  analyze  quantitatively  the  effect  of  flame  collisions  and  to  determine  the  actual  increase  of  S  /  in  the 
resulting  cusps,  consider  an  idealized  model  of  the  cusps  observed  in  Fig.  9.  In  particular,  consider  two  symmetrically 
located  planar  flame  sheets  approaching  each  other  with  the  inclination  angle  to  the  mid-plane  a.  At  the  point  of 
collision,  the  flame  sheets  merge  and  form  a  cusp.  Its  properties,  such  as  its  structure,  speed  of  propagation,  and 
effective  burning  velocity,  are  determined  in  this  configuration  only  by  two  parameters,  S  l  and  a. 

Fig.  10  illustrates  the  flame  structure  formed  in  this  situation.  It  shows  the  distribution  of  the  fuel  mass  fraction  in 
two  simulations  performed  for  a  =  1°  and  a  -  4°  with  Athena-RFX  using  the  same  physical  model  as  in  simulations 
S 1  -  S3.  The  domain  has  the  resolution  Ax  =  S/J32  and  zero-order  extrapolation  boundary  conditions  on  all  sides.  At 
t  =  0,  two  intersecting  flame  sheets  were  initialized  with  the  exact  structure  of  the  planar  laminar  flame,  and  uniform 
pressure  Pq  and  zero  velocities  were  assumed.  After  the  initial  transient  stage,  the  flame  develops  the  structure  shown 
in  Fig.  10.  The  structure  of  the  cusp  itself  does  not  change  with  time,  provided  that  the  planar  flame  sheets  extend 
sufficiently  far  from  it.  Therefore,  Fig.  10  can  be  viewed  as  a  steady-state  solution. 

Fig.  10  shows  that  flame  collision  at  very  low  inclination  angles  results  in  the  formation  of  a  highly  elongated 
structure.  It  is  formed  by  rapid  heating  and,  subsequently,  ignition  of  a  larger  amount  of  fuel  than  in  the  planar 
laminar  flame  due  to  the  focusing  of  the  thermal  flux  from  two  approaching  flame  sheets.  Consequently,  its  extent  is 
determined  by  the  size  of  the  region  in  which  the  preheat  zones  of  approaching  flames  overlap  and,  thus,  create  the 
necessary  enhanced  thermal  flux.  Since  the  width  of  the  preheat  zone  is  « ()/.,  very  small  values  of  a  are  required  for 
the  cusp  to  be  substantially  broader  than  8i .  This  can  further  be  seen  in  the  distributions  of  the  fuel  mass  fraction, 
Y,  temperature,  7  ,  and  reaction  rate,  Y,  along  the  symmetry  axis  of  the  cusp  shown  in  Fig.  11.  The  structure  of  the 
reaction  zone  in  the  cusp  very  quickly  approaches  that  of  the  laminar  flame  with  two  becoming  very  close  to  each 
other  already  at  a  ~  4°.  At  the  same  time,  the  preheat  zone  of  the  cusp  remains  substantially  broadened  for  much 
larger  values  of  a. 

There  are  three  characteristic  speeds  in  this  problem.  Sufficiently  far  from  the  tip  of  the  cusp,  the  flame  locally 
propagates  normal  to  its  surface  with  the  laminar  flame  speed  S / ,  i.e.,  there  Si  —  S l-  The  cusp  itself  moves  into  the 
fuel  with  the  phase  velocity  Dc,  which  results  simply  from  the  two  inclined  planar  surfaces  colliding  with  each  other. 
This  speed  can  be  determined  based  on  purely  geometrical  considerations  [25], 

Dc  =  —.  (20) 

sin  cr 

This  cusp  propagation  speed,  normalized  by  S l,  is  shown  as  a  solid  line  in  Fig.  12a  along  with  Dc/S l  determined 
as  the  velocity  of  the  leftmost  point  of  the  Y  -  0.5  isosurface  in  simulations  for  four  values  of  a.  Eq.  (20)  is  within 
<  1%  accuracy  of  the  computed  values.  Dc  approaches  the  laminar  flame  speed  as  a  — »  90°,  i.e.,  when  the  two  flame 
sheets  cease  to  advance  toward  each  other.  At  the  opposite  limit  of  small  a,  Dc  can  be  quite  large  and,  in  principle, 
infinite  when  a  — »  0°.  In  particular,  at  a  ~  0.5°,  l),  becomes  larger  than  the  sound  speed  in  cold  fuel  indicated  with 
the  horizontal  dashed  line  in  Fig.  12a.  Note  also  that  at  a  =  1°,  the  value  of  I),  »  57. S’/,  is  very  close  to  the  high 
propagation  velocity  &55S  l  of  cusps  in  regions  B  -  C  in  Fig.  9  estimated  in  §  5.2. 

The  third  characteristic  speed,  which  is  of  most  importance  to  us,  is  the  local  burning  speed  in  the  cusp,  S  /.  The 
broadened  reaction  zone  at  low  flame-inclination  angles  (Figs.  10  and  1  lb)  causes  more  fuel  to  be  consumed  per  unit 
flame  surface  area  than  in  the  planar  laminar  flame,  and,  therefore,  it  leads  to  a  larger  local  flame  speed  S  /  >  S  l-  Since 
the  reaction  zone  width  is  the  largest  at  the  tip  of  the  cusp,  S  /  has  its  maximum  there  and  it  gradually  decreases  to  its 
laminar  value  S  l  in  the  planar  regions  of  the  flame  sufficiently  far  from  the  cusp.  The  maximum  value  of  S  /  at  the 
cusp  tip  can  be  found  as 


SI  =  max(S /)  =  —  f pYdx  =  f  p2Y  exp  (  — — )  dx.  (21) 

po  J  Po  J  \  RT I 

Here  eq.  (7)  was  used  for  Y ,  and  the  integral  is  taken  along  the  symmetry  axis  of  the  cusp  shown  in  Fig.  10.  Fig.  12b 
shows  the  computed  values  of  S  J  normalized  by  S l  for  the  same  four  inclination  angles  as  in  Fig.  12a. 

Two  key  conclusions  emerge  from  Fig.  12.  First,  S*  is  substantially  lower  than  Dc,  being  less  by  more  than  an 
order  of  magnitude  at  the  values  of  a  considered.  Second  and  most  important,  S  *  can  indeed  be  substantially  higher 
than  S l  at  small  inclination  angles. 
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5.4.  Connection  between  the  increase  ofS  /  in  cusps  and  S  r 

We  are  interested,  however,  not  just  in  S  /  but,  rather,  in  the  burning  speed  of  an  extended  flame  region,  as  this  is  a 
direct  equivalent  of  the  turbulent  flame  speed.  The  total  burning  speed  of  the  flame  configuration  shown  in  Fig.  10  is 


o  C  —  - , 

P0 

where  Mr  is  the  total  mass  of  reactants  converted  into  product  per  unit  time.  If  the  flame  were  to  propagate  everywhere 
with  the  speed  S l,  then  Sc  =  AS  l,  where  A  is  the  surface  area  of  the  flame  in  such  a  configuration  based  on  the 
isosurface  of  Y  =  0.15  representing  the  peak  reaction  rate.  Otherwise, 


I  = 


Sc 

ASl 


(23) 


Here  I  again  is  the  stretch  factor,  which  is  equivalent  to  the  definition  used  in  eq.  (12)  in  the  context  of  the  turbulent 
flame  speed.  Substituting  the  definition  of  S  t  given  by  eq.  (10)  into  the  definition  of  I  given  by  eq.  (12),  it  can  be  seen 
that  indeed  m;r/Po  —  IAS l  =  S c. 

Since  S  /  monotonically  decreases  away  from  the  cusp,  both  S ,  and  I  depend  on  the  size  of  the  flame  region  being 
considered.  This  size  can  be  represented  by  the  length  lc,  illustrated  in  Fig.  10,  which  is  defined  as  the  distance  from 
the  leftmost  point  of  the  isosurface  of  Y  —  0.15  in  the  direction  of  cusp  propagation.  Considering  different  values  of  lc 
is  equivalent  to  varying  the  fraction  of  the  total  surface  area  of  the  flame  in  Fig.  10  represented  by  the  planar  section 
in  which  Si  ~  S l.  Consequently,  this  allows  one  to  vary  the  relative  contribution  to  S c,  and  thus  to  /,  of  the  cusp  in 
which  S i  >  Sl. 

The  stretch  factor,  calculated  using  eqs.  (22)-(23),  is  shown  in  Fig.  13b  for  four  values  of  a  as  a  function  of  Ic/6l. 
At  smaller  lc  close  to  the  width  of  the  reaction  zone  for  a  given  a  (Fig.  1  lb),  /  essentially  reflects  only  the  values  of 
S  i  in  the  immediate  vicinity  of  the  tip  of  the  cusp.  On  the  other  hand,  at  larger  lc,  the  surface  area  of  the  planar  flame 
region,  in  which  Si  «  S l,  is  larger  and,  thus,  such  region  represents  a  greater  fraction  of  the  total  surface  area  of  the 
flame  configuration  shown  in  Fig.  10.  As  a  result,  the  contribution  of  the  planar  flame  starts  to  dominate  that  of  the 
cusp,  which  causes  /  — »  1  as  seen  in  Fig.  13b.  At  lower  inclination  angles,  however,  the  region  in  which  S  /  >  S  / 
extends  progressively  further  from  the  cusp  tip.  In  particular,  at  a  =  1°,  /  is  substantially  greater  than  unity  even  at 
Ic  »  6l.  At  the  same  time,  already  at  a  =  4°,  I  drops  to  <  1.04  at  lc  «  56 l  which  shows  that  S /  recovers  its  laminar 
value  within  a  few  laminar  flame  widths  from  the  cusp  tip. 

These  results  show  that  the  fomiation  of  a  cusp  due  to  the  collision  of  planar  flame  sheets  indeed  produces  values 
of  /  substantially  larger  than  unity  and  comparable  to  those  observed  in  the  simulations  presented  here.  As  a  result 
of  a  higher  local  flame  speed  in  a  cusp,  its  contribution  to  the  global  turbulent  burning  speed  is  disproportionately 
large  compared  to  the  fraction  of  the  instantaneous  flame  surface  area  in  it.  A  more  complex  flame  configuration, 
consisting  of  multiple  planar  flame  sheets  that  approach  each  other  and  collide  forming  multiple  cusps,  can  then  be 
viewed  as  a  simple  model  of  the  turbulent  flame.  If  the  flame  surface  density  in  this  system  is  increased,  flame  sheets 
will  get  closer  to  each  other  causing  flame  collisions  to  become  more  frequent  and  the  resulting  cusps  to  become 
more  numerous.  Consequently,  the  fraction  of  the  flame  surface  area  contained  in  cusps  will  increase  and  I  will  grow 
rapidly  and  significantly,  which  is  demonstrated  by  the  increase  in  /  with  decreasing  Ic. 

This  simplified  model  does  not  take  into  account  the  full  complexity  of  an  actual  turbulent  flame.  Such  flame  does 
not  consist  of  perfectly  planar  flame  sheets  that  merge,  but  instead  it  is  constantly  wrinkled  and  folded  on  a  variety 
of  scales.  The  analysis  presented  in  §  5.3  shows  that  in  a  curved  flame,  S  /,  unlike  Dc,  does  not  gradually  rise  with 
increasing  curvature.  Instead,  S  /  remains  very  close  to  its  laminar  value  until  the  curvature  becomes  very  large,  i.e., 
>  I  / 6 1_,  which  is  consistent  with  the  behavior  of  an  ideal  collapsing  spherical  or  cylindrical  flame  discussed  in  §  5.2. 
Turbulence  can  create  regions  of  such  large  curvature  either  by  folding  the  flame,  thereby,  gradually  increasing  the 
curvature  at  a  specific  point  of  the  flame  surface  until  it  becomes  >  I  /()/,,  or  by  bringing  two  flame  sheets  together 
until  they  merge.  Both  of  these  situations,  however,  are,  in  fact,  flame  collisions.  When  they  occur,  the  cusp  forms 
and  only  then  the  local  flame  speed  in  it  increases. 

Based  on  these  considerations,  the  following  picture  emerges.  At  low  turbulent  intensities,  the  flame  surface  is 
wrinkled  (or  folded)  primarily  on  larger  scales  with  the  curvature  radius  at  each  point  »  6/ .  Consequently,  the  flame 
propagates  locally  with  its  laminar  speed,  and  the  total  St  °c  At.  Individual  cusps  can  form  even  at  low  turbulent 
intensities  as  a  result  of  the  collision  of  individual  flame  sheets.  At  low  values  of  2,  however,  cusps  are  infrequent 
and  they  represent  only  a  small  fraction  of  the  total  flame  surface  resulting  in  I  «  1 .  This  is  analogous  to  the  situation 
considered  above  when  lc  is  large.  As  the  turbulent  intensity  increases,  turbulence  is  not  able  to  disrupt  the  internal 
structure  of  the  flamelets  causing  flamelets  to  retain  locally  their  laminar  structure  and  velocity.  Instead,  turbulent 
motions  are  much  more  efficient  at  folding  the  flame  with  increasingly  greater  curvature,  creating  a  tightly  packed 
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configuration  with  large  £.  This  increases  the  frequency  of  flame  collisions  and  results  in  a  much  higher  rate  of  cusp 
formation.  Consequently,  a  progressively  larger  fraction  of  Aj  is  contained  in  cusps,  which,  by  analogy  with  lower 
values  of  lc,  leads  to  I  >  1. 

This  picture  reconciles  the  seeming  contradiction  discussed  above  between  the  need  for  the  enhanced  local  flame 
speed,  as  suggested  by  values  of  /  >  1,  and  the  absence  of  such  enhancement  globally  throughout  the  flame  by  the 
turbulent  diffusive  transport,  as  evidenced  by  the  absence  of  the  flame  broadening.  The  flame  speed  can  be  significantly 
enhanced  locally  in  cusps  due  to  thermal  flux  focusing  and  not  because  the  flame  itself  is  disrupted  and  broadened  by 
turbulence. 

This  simple  model  of  cusp  formation  demonstrates  the  mechanism  through  which  5  j  can  produce  an  exaggerated 
response  to  the  increase  in  At.  At  the  same  time,  in  an  actual  turbulent  flame  cusps  will  form  in  a  multitude  of 
configurations.  For  instance,  in  addition  to  the  purely  two-dimensional  situation  considered  here,  there  will  also  be 
3D  flame  collisions  which  will  result  in  even  higher  local  flame  speeds  in  the  cusps  and,  thus,  will  have  an  even 
larger  effect  on  the  magnitude  of  I.  Therefore,  in  order  to  predict  a  particular  value  of  I  which  can  be  expected  in  the 
turbulent  flow  of  a  specific  intensity,  it  is  necessary  to  understand  the  types  of  cusps  which  can  form,  the  local  flame 
speed  of  each  type,  its  probability  of  formation,  and,  thus,  the  contribution  of  each  type  to  the  exaggerated  response 
ofSr.  Such  detailed  analysis  is  the  subject  for  future  studies. 


5.5.  Criterion  for  onset  of  the  cusp-dominated  regime  of  flame  evolution 

The  model  of  cusp  formation  can  be  used  to  determine  the  critical  turbulent  intensity  above  which  the  evolution 
of  a  turbulent  flame  can  be  expected  to  be  dominated  by  cusps,  thereby,  leading  to  values  of  I  substantially  larger 
than  unity.  Consider  a  section  of  an  idealized  curved  flame  front  containing  a  cusp  and  perturbed  with  a  wavelength 
Ac  and  an  amplitude  f,  schematically  shown  in  Fig.  14.  This  type  of  structure  was  considered  by  Zel’dovich  [25]  as 
a  model  of  the  flame  formed  under  the  action  of  the  Landau-Darrieus  instability  in  order  to  analyze  the  stabilizing 
effect  of  cusps.  Such  idealized  flame,  however,  can  develop  under  the  action  of  any  destabilizing  process,  e.g.,  the 
Rayleigh-Taylor  instability  or,  in  our  case,  turbulence.  In  particular.  Fig.  14  can  represent  the  flame  surface  deformed 
by  two  adjacent  counter-rotating  vortices  of  size  Ac. 

As  discussed  in  §  5.3  (Fig.  12),  the  speed  of  cusp  propagation,  Dc,  (Fig.  14)  is  much  larger  than  both  S  l  and  the 
maximum  local  flame  speed  in  the  cusp,  S  ,* .  Moreover,  unlike  Sf  Dc  »  S  /  even  at  large  flame  inclination  angles, 
and,  thus,  it  increases  gradually  with  curvature.  Therefore,  Dc  is  the  primary  factor  responsible  for  decreasing  the 
amplitude  of  the  cusp  and  smoothing  the  flame  surface.  In  particular,  the  rate  of  decrease  of  the  cusp  amplitude  is  [25] 


(24) 


The  angle  a  and,  thus,  ( dlfldt )_,  will  change  with  the  cusp  amplitude.  In  order  to  relate  a  and  lc,  a  parabolic  shape  of 
the  flame  was  assumed  in  [25],  and  eq.  (24)  can  then  be  rewritten  as 


(25) 


The  net  rate  of  change  of  lc  is  determined  by  the  balance  of  the  stabilizing  process  described  by  eq.  (25)  and  a 
destabilizing  one,  which  gives 


die 

dt 


=  x¥  -8S 


(26) 


In  the  turbulent  flow,  flame  perturbations  can  be  assumed  to  grow  linearly  in  time  as  they  are  being  stretched  by  the 
turbulent  speed  U characteristic  of  the  scale  Ac.  This  gives  the  growth  rate  T  =  UA.  By  setting  dlfldt  -  0,  the 
limiting  value  of  lc(Ac)  can  be  determined.  It  was  shown  by  Khokhlov  [7]  in  the  context  of  Rayleigh-Taylor-unstable 
flames  that  nonlinear  flame  stabilization  due  to  cusp  propagation  can  be  expected  to  fail  once  lc  >  Ac.  Therefore,  by 
setting  lc  =  Ac,  'F  =  U,\,  and  dlfldt  —  0  in  eq.  (26),  he  found  the  critical  value  of  the  turbulent  velocity  on  a  given 
scale  Ac  [7] 

U,  =  8  SL.  (27) 


This  is  the  value  of  U,\  that  is  needed  to  overcome  the  stabilizing  effect  of  cusps  and,  thereby,  to  deform  the  flame  on 
scale  Ac.  This  result  was  used  to  demonstrate  that,  in  order  for  the  flame  to  be  unstable  at  a  wavelength  Ac  under  the 
action  of  turbulent  motions,  the  turbulent  speed  on  that  scale  must  be  substantially  larger  than  S  l  [7]. 

This  shows  that  given  a  specific  turbulent  cascade,  the  flame  will  be  smooth  on  scales  smaller  than  some  critical 
wavelength  A*  on  which  the  condition  in  eq.  (27)  is  satisfied.  It  was  shown  above  that  in  order  for  the  cusps  to  have 
a  pronounced  effect  on  the  flame  evolution,  the  flame  must  be  folded  on  scales  comparable  to  the  full  flame  width  If, 
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which  in  our  case  is  ~  26  F.  Therefore,  by  setting  A*  =  lF  and  using  eq.  (27),  we  find  that  the  flame  will  be  curved 
on  scales  ~lF  when  turbulent  velocity  on  that  scale  becomes  ~  85/..  In  a  turbulent  flow  with  the  Kolmogorov  energy 
spectrum,  the  corresponding  critical  value  of  the  integral  turbulent  velocity  then  is 


u;  ~  85 


(28) 


This  gives  the  critical  values  of  the  Karlovitz  number,  Ka* ,  and  the  Damkohler  number,  Da*,  written  using  their 
traditional  definitions  [10] 


Ka * 
Da* 


—  -  l—\l2 

tn  VTc/ 
h  _  IS  l  _ 
If  IfUi 


=  8 3/2  »  20, 

1(1)* 

8 \lF) 


(29) 


Here  tF  =  IF/S l  is  the  characteristic  flame  time,  tn  is  the  Kolmogorov  time,  tF  =  I/Ui  is  the  characteristic  turbulent 
time  on  the  integral  scale,  and  Lq  =  I(S l/Ui)2  is  the  Gibson  scale. 

The  range  of  the  regimes,  in  which  I  is  expected  to  be  larger  than  unity  according  to  the  eq.  (28),  is  shown  on  a 
traditional  combustion  regime  diagram  [10]  in  Fig.  15  as  the  orange  region.  The  filled  red  square  corresponds  to  the 
simulations  presented  here.  We  also  determined  the  time-averaged  value  of  /  in  a  simulation  with  approximately  twice 
lower  turbulent  intensity  which,  however,  was  still  above  the  Ka  «  20  line.  In  this  calculation,  which  is  represented 
with  the  open  red  square  in  Fig.  15,  we  found  7 0.15  ~  1.16.  Details  of  this  calculation  will  be  presented  in  a  separate 
paper. 

Despite  the  simplicity  of  the  cusp  model  used  here,  eqs.  (28)-(29)  provide  a  rather  accurate  criterion  for  the  onset 
of  the  cusp-dominated  regime  of  the  turbulent  flame  evolution.  In  particular,  the  decrease  in  the  value  of  7o.i5  between 
the  two  turbulent  intensities  considered  in  Fig.  15  suggests  that  below  the  Ka  ~  20  line,  I  deviates  from  unity  at  most 
by  a  few  percent.  Consequently,  this  line  can  be  viewed  as  an  approximate  upper  range  of  validity  of  Damkohler’s 
concept. 


5.6.  Nonlinear  regime  of  turbulent  flame  evolution 


Finally,  how  will  the  distribution  of  5  7-/5 /.  as  a  function  of  Aq^/L2  shown  in  Fig.  5b  vary  with  the  change  in 
turbulent  intensity? 

Based  on  the  model  discussed  in  §  5. 3-5.4,  the  total  turbulent  flame  speed  can  be  viewed  as  consisting  of  two 
components.  The  first  is  the  contribution  of  the  smoothly  curved  flame  regions  with  the  curvature  radius  >6l  in  which 
5/  =  5  l-  The  second  is  the  effect  of  cusps.  Since  5/  >  5/.  in  cusps,  their  contribution  depends  nonlinearly  on  the 
flame  surface  area  contained  in  them,  as  evidenced  by  Figs.  13. 6  Therefore,  the  total  fuel-consumption  speed  of  the 
turbulent  flame  can  be  written  as  ivr/po  =  (1  -  fc)S  lAj  +  fcS iflF' CAy),  where  T'lAj)  is  some  nonlinear  function,  and 
fc  is  the  fraction  of  the  flame  surface  area  contained  in  cusps.  This  can  be  rewritten  in  terms  of  5  7-  defined  in  eq.  (10) 
as 


5  T 
5l 


=  a -/,)§  + 


(30) 


where  f  and  T  are  nonlinear  functions. 

The  balance  between  the  linear  and  nonlinear  terms  in  eq.  (30)  is  controlled  by  fc  and,  therefore,  by  the  rate 
of  cusp  creation,  ( d(fcAF)/dt)+ .  The  following  qualitative  physical  model  can  be  used  to  estimate  this  rate  (see 
also  [7]).  Consider  a  structure  consisting  of  planar  flame  sheets,  as  was  discussed  in  §  5.4.  If  the  flame  surface 
density  of  such  configuration  is  Xy ,  then  the  average  flame  separation  is  1  /Xy  .  Such  flame  sheets  will  move  toward 
each  other  with  the  speed  5  l  +  U t,  where  Uj  is  some  characteristic  turbulent  velocity  responsible  for  the  advective 
transport  of  the  flame.  Since  Ui  is  the  largest  velocity  of  coherent  turbulent  motions,  then  UF  can  be  approximated 
with  Ui.  Consequently,  the  flame  sheets  will  merge  and  annihilate  in  the  time  tc  ~  1  /CZt(Sf  +  Ui)).  The  quantity 
(d(fcAF)/dt)+  will  depend  on  the  total  surface  area  of  the  flame  sheets  times  the  frequency  of  their  collisions,  i.e., 
( d(fcAF)/dt)+  oc  Aj/tc  oc  Aj'LjiS l  +  Uf).  While,  in  principle,  the  change  in  X7-  does  not  imply  a  similar  change  in 
Aj,  and  vice  versa,  it  was  shown  in  §  4.3  that  AF  and  Xy-  are  well  correlated.  Therefore,  Aj  can  be  viewed  as  simply 
proportional  to  X7  and,  thus,  finally 


(d(fd^T>)+  (Sl  +  UipT  (Sl  +  U)A2. 


(31) 


6Note  that  the  flame  surface  area  of  the  configuration  shown  in  Fig.  10  does  not  depend  on  a  and  is  a  function  only  of  lc. 
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Note  that  this  expression  is  very  similar  to  the  destruction  term  often  used  in  the  balance  equation  of  the  flame  surface 
density  (see  [28]  for  a  review  of  a  number  of  such  models).  This  is  also  analogous  to  the  rate  of  collisions  of  gas 
molecules  with  Zy  playing  the  role  of  number  density  and  the  turbulent  intensity  playing  the  role  of  temperature  in 
determining  the  speed  with  which  constituents  approach  each  other. 

In  §  5.4  we  noted  that  in  the  actual  turbulent  flame,  there  will  exist  other  types  of  flames  collisions  besides  the  ones 
that  can  be  represented  with  planar  flame  sheets.  Such  configurations  will  result  in  an  even  faster  rate  of  cusp  creation 
that  would  likely  be  proportional  to  the  higher  powers  of  Zy  and,  possibly,  U /.  The  probability  of  the  formation  of 
such  flame  collisions,  however,  will  rapidly  decrease  with  the  increase  in  their  complexity  [3,  7]. 

Eq.  31  shows  that  the  flame  surface  area  contained  in  cusps  in  the  turbulent  flame  grows  faster  than  the  total  Ay. 
Consequently,  with  increasing  Ay  and  Zy,  the  balance  in  eq.  (30)  shifts  rapidly  toward  the  nonlinear  term  associated 
with  cusps. 

This  picture  then  shows  how  the  distribution  of  Sy/Sy  as  a  function  of  A/L 2,  given  in  Fig.  5b,  changes  with 
turbulent  intensity.  At  a  specific  [//,  a  variety  of  different  flame  configurations  with  different  values  of  Ay  and  Zy 
are  realized  in  the  course  of  flame  evolution.  Larger  values  of  Ay  lead  to  larger  Zy,  which  increases  the  frequency 
of  occurrence  of  cusps,  i.e.,  fc.  As  a  result,  the  nonlinear  term  in  eq.  (30),  which  rapidly  grows  with  Ay,  begins  to 
dominate,  and  this  causes  a  progressively  more  exaggerated  response  of  S  y  to  the  increase  in  Ay.  Moreover,  since  a 
flame  with  a  given  surface  area  can  exist  in  a  variety  of  configurations  with  different  numbers  and  types  of  cusps,  this 
creates  a  scatter  in  the  distribution  of  Sy/Sy  vs.  Ay/L2. 

At  lower  turbulent  speeds,  the  distribution  of  Sy/Sy  vs.  Ay/L 2  should  collapse  onto  the  Sy/Sy  =  Ay/L2  line  since, 
in  the  limit  of  complete  absence  of  cusps,  there  is  a  unique  correspondence  between  these  two  quantities.  Eq.  (30) 
suggests  that  such  behavior  would  indeed  take  place.  As  V /  becomes  smaller,  the  flame  will  be  less  convolved,  causing 
the  center  of  mass  of  the  distribution  to  shift  toward  smaller  values  of  Ay/L 2  and,  thus,  Sy/Sy.  Since  the  surface  area 
contained  in  cusps  depends  not  only  on  Ay  but  also  on  U/  (eq.  (31)),  smaller  [//  will  cause  the  cusps  to  be  less  prevalent 
at  a  fixed  Ay/L2.  As  a  result,  the  nonlinear  term  in  eq.  (30),  will  be  smaller  and  the  overall  distribution  Sy/Sy  vs. 
Ay/L2  will  be  closer  to  the  Sy/Sy  =  Ay/L 2  line.  Moreover,  since  all  potential  statistical  variations  in  the  distribution 
are  associated  with  cusps  and  are,  thus,  represented  by  the  nonlinear  term  TiA/L2),  the  decrease  in  the  latter  will  also 
lead  to  the  decrease  in  the  distribution  scatter. 

At  larger  Uy  the  behavior  will  be  opposite.  The  center  of  mass  of  the  distribution  will  shift  to  larger  values  of 
Ay/L 2  and  Sy/Sy.  The  nonlinear  term  in  eq.  (30)  will  become  progressively  more  dominant,  causing  larger  values  of 
Sy/Sy  at  a  given  Ay  /L2  and,  thus,  enhancing  the  exaggerated  response  of  5  y . 

As  the  turbulent  intensity  increases,  the  overall  distribution  of  Sy/Sy  vs.  Ay/L2  will  not  only  deviate  more 
strongly  from  the  Sy/Sy  -  Ay  /L2  limit,  but  it  will  also  become  more  pronouncedly  curved  upward.  This  is  the  result 
of  the  fact  that 


d2Sy(Ay ) 
dA 


(32) 


and  the  magnitude  of  d2S  y(Ay)/dAj.  at  every  point  increases  as  U\  becomes  larger.  Here  Sy(Ay)  represents  the 
time-averaged  value  of  S  y  in  the  limit  of  large  statistics. 

To  show  that  inequality  (32)  holds,  consider  first  the  limiting  value  of  Zy.  The  fact  that  the  flame  is  not  an 
infinitely  thin  surface  means  that  such  limit  does  exist.  Once  Zy  becomes  large  enough  so  that  the  individual  flame 
sheets  come  into  contact,  Zy  cannot  increase  beyond  that.  This  maximum  value  Z yjnax  can  be  estimated  assuming 
that  the  minimum  flame  separation  is  equal  to  the  full  flamelet  width  /y,  which  in  our  case  is  ~  26 y  (cf.  Fig.  7  in  [9]). 
Then7 

S'T.max  —  Eq.  1 5. max  —  ~  1-5  mm  .  (33) 

2  by 

It  was  discussed  in  [9],  based  on  the  comparison  of  the  results  of  that  work  with  the  results  of  Aspden  et  al.  [29],  that 
the  turbulent  flame  width,  6y,  is  primarily  determined  by  the  integral  scale,  /,  (and,  therefore,  by  the  driving  scale,  L). 
A  similar  conclusion  was  also  reached  by  Peters  [10]  (cf.  eq.  (2.175)  therein).  If  L  is  held  fixed,  then  /  will  also  remain 
constant  as  the  turbulent  intensity  increases.  Consequently,  6y  will  not  change  with  increase  in  U /.  Substituting  the 
value  of  6y  obtained  in  simulation  S3  (see  Table  3)  into  eq.  (18),  and  setting  Zy  =  Z y.max,  the  corresponding  maximum 
value  of  the  flame  surface  area  is 

At, max  _  ^0.15 .max  a  \ 

~U~  =  L2  *  (  ) 

Part  of  the  increase  in  Ay  is  typically  associated  with  the  increase  in  6y  and  not  only  Zy.  Therefore,  actual  values  of 
Ay, max/ L2  can  be  somewhat  larger  than  the  one  given  above.  Thus  eq.  (34)  should  be  viewed  as  a  more  approximate 
limit  than  eq.  (33). 


7Note  that  Eymai  is  different  from  Z.max  typically  used  in  combustion  research  and  discussed  in  §  4.1. 
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The  situation  when  all  of  the  flame  surface  is  in  contact  with  itself  corresponds  to  the  infinite  local  flame  speed 
and,  thus,  infinite  5 t-  Therefore,  as  Ay / L 2  — »  Ax.max/L2,  both  St/S l  —■ >  00  and  clS T  /d  Ay  — >  00.  On  the  other  hand, 
at  small  values  of  At/L2,  the  slope  of  S  t  is  finite  and,  in  fact,  as  At /L2  — >  1 ,  dS  t /dAT  — >  S 1/L2.  This  shows  that 
the  slope  of  the  distribution  St/Sl  cannot  be  constant  and  instead  dS  7  / d/\  7  has  to  be  a  monotonically  increasing 
function,  which  proves  the  inequality  (32).  Since  this  inequality  is  the  consequence  of  the  creation  of  cusps  in  the 
turbulent  flame,  the  magnitude  of  d2S  TiA^/dA^  will  increase  with  increase  in  the  frequency  of  creation  of  cusps 
and,  thus,  with  increase  in  both  At  and  U\.  Therefore,  as  turbulent  intensity  becomes  larger,  the  curvature  of  the 
distribution  will  become  pronounced  at  progressively  lower  values  of  At. 

6.  Conclusions 

This  work  continued  the  analysis  of  the  set  of  three  numerical  simulations  first  presented  in  [9],  These  calculations 
model  the  interaction  of  the  premixed  flame  with  high-speed,  subsonic,  homogeneous,  isotropic  turbulence  in  an 
unconfined  system,  i.e.,  in  the  absence  of  walls  and  boundaries.  The  turbulent  r.m.s.  velocity,  Urms,  is  ~  35  times  larger 
than  the  laminar  flame  speed,  S  l-  The  resulting  Damkohler  number  based  on  the  turbulent  integral  scales  is  Da  =  0.05. 
It  was  demonstrated  in  [9]  that  this  system  represents  turbulent  combustion  in  the  thin  reaction  zone  regime.  Even  in 
the  presence  of  such  intense  turbulence,  the  turbulent  flame  consists  of  highly  convolved  flamelets  with  the  reaction- 
zone  structure  virtually  identical  to  that  of  a  planar  laminar  flame  and  with  the  preheat  zone  broadened  by  a  factor 
*2. 

The  fact  that  turbulence  is  unable  to  penetrate  and  disrupt  the  internal  flame  structure  showed  that  diffusive  pro¬ 
cesses  are  not  enhanced  on  small  scales  by  turbulent  transport,  and  flamelets  propagate  locally  with  the  speed  of  the 
planar  laminar  flame  [9].  This  raised  the  following  question:  Can  the  magnitude  of  the  turbulent  flame  speed  in  the 
presence  of  high- intensity  turbulence  be  fully  accounted  for  by  the  increase  in  the  flame  surface  area,  as  was  originally 
suggested  by  Damkohler  [1]?  Here  we  summarize  the  main  findings  of  our  study  of  this  issue. 

Analysis  of  the  area  and  density  of  the  fuel  mass-fraction  isosurfaces,  A{Y)  and  X(K ),  showed  that  the  flamelet, 
folded  inside  the  flame  brush  in  the  presence  of  high-speed  turbulence,  cannot  be  viewed  as  a  thin  uniform  structure, 
in  which  all  isosurfaces  are  parallel  to  each  other.  Different  regions  of  the  flamelet  have  quite  different  response  to  the 
action  of  turbulence.  In  the  higher-resolution  calculations  S2  and  S3,  both  A  and  X,  on  average,  increase  monotonically 
through  the  flamelet  with  decreasing  temperature,  which  is  manifested  in  the  distinctive  inverted-S  shape  of  their  time- 
averaged  distributions.  As  a  result,  isosurfaces  of  higher  fuel  mass  fractions  are  folded  by  turbulence  on  progressively 
smaller  scales.  This  causes  the  substantially  finer  wrinkling  of  the  flame  surface  on  the  fuel  side  than  on  the  product 
side  observed  in  [9]. 

Distributions  of  A(Y)  and  X(K)  showed  that  in  the  presence  of  the  high-speed  turbulence  considered  here,  the 
definition  of  the  flame  surface  area.  At,  must  be  revisited  before  the  balance  between  A  j  and  the  turbulent  flame 
speed,  St,  can  be  considered.  In  particular,  it  must  be  determined  which  value  of  Y  characterizes  the  evolution  and 
global  properties  of  the  turbulent  flame  most  fully  and  accurately.  In  order  to  answer  this  question,  we  analyzed  the 
correlation  of  A(Y)/L2  with  quantities  that  characterize  both  the  energetics  and  the  global  structure  of  the  turbulent 
flame,  namely  its  normalized  speed,  St/Sl,  and  width,  6t/L.  We  also  analyzed  the  correlation  of  X(  Y )  with  St/Sl- 
This  analysis  demonstrated  that  global  properties  of  the  turbulent  flame  are  best  represented  by  the  structure  of  the 
region  of  peak  reaction  rate.  For  the  reaction-diffusion  model  used  in  this  work,  this  corresponds  to  Y  «  0.15. 
Therefore,  the  isosurface  of  this  value  of  Y  must  be  viewed  as  the  flame  surface  and  At  =  A0.15. 

Larger  values  of  At  are  associated  with  larger  flame  surface  density,  Xy .  In  other  words,  Ay-  grows  primarily 
as  a  result  of  the  much  tighter  folding  and  packing  of  the  flame,  rather  than  only  due  to  the  increase  in  the  overall 
width  of  the  flame  brush.  Given  the  absence  of  any  broadening  of  the  reaction  zone  observed  in  [9],  this  shows  that 
high-intensity  turbulence  is  much  more  efficient  at  tightly  packing  the  flamelets  inside  the  flame  brush  rather  than  at 
disrupting  and  broadening  their  internal  structure  (also  see  [10]). 

Consideration  of  the  stretch  factor,  I,  calculated  at  Y  =  0.15,  and,  more  generally,  of  the  distribution  of  St/Sl  vs. 
A0.15/L2  led  to  the  following  key  conclusion  of  this  work:  In  the  presence  of  high-speed  turbulence,  the  magnitude 
of  the  turbulent  flame  speed  cannot  be  attributed  purely  to  the  increase  in  the  flame  surface  area.  In  particular,  in 
the  highest-resolution  simulation  S3  at  all  times  /0.15  >  1  and,  thus,  St/Sl  >  A0.15/L2.  The  deviation  of  /0.15  from 
unity  is,  on  average,  «30%  and  occasionally  reaches  as  high  as  50%.  Since  the  local  flame  speed  is  not  increased  by 
turbulence  [9],  this  shows  that  Damkohler’s  concept  breaks  down  for  sufficiently  high-intensity  turbulence,  even  for 
the  flows  characterized  by  Le  —  1 . 

The  deviation  of  S  t  /S  l  from  At/L1  becomes  larger  as  Ay-  and  Xy  increase.  In  other  words,  an  increase  in  the  flame 
surface  area,  and  the  associated  tighter  packing  of  the  flame,  result  not  just  in  a  larger  turbulent  flame  speed  but,  most 
importantly,  in  the  progressively  more  exaggerated  response  of  St  to  the  increase  in  Ay.  This  causes  substantially 
accelerated  burning. 
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Tightly  packed  flame  configurations,  produced  by  high-speed  turbulence,  result  in  frequent  flame  collisions,  which 
lead  to  the  formation  of  regions  of  high  flame  curvature  >  1  / 6/ ,  or  “cusps.”  This  results  in  significant  focusing  of 
the  thermal  flux  over  an  extended  region  of  the  flame  surface  which  increases  the  local  flame  speed  in  the  cusp,  S  /, 
over  its  laminar  value,  S  l-  Due  to  the  large  values  of  S  /,  the  contribution  of  cusps  to  the  total  S  y  is  disproportionately 
large  compared  to  the  flame  surface  area  in  them.  This  provides  a  natural  mechanism  for  the  formation  of  the  exag¬ 
gerated  response  of  S  t ■  The  increase  of  S  /  in  cusps  is  inherently  local,  and  it  does  not  require  flame  broadening  and 
acceleration  by  turbulent  transport,  in  agreement  with  the  results  of  [9]. 

Our  results  suggest  that  there  exist  two  distinct  regimes  of  flame  evolution.  At  low  turbulent  intensities,  the 
turbulent  flame  evolves  in  the  linear  regime.  Here  the  role  of  cusps,  and  thus  of  the  nonlinear  term  in  eq.  (30),  is 
negligible  and  St  °c  Aj.  High  turbulent  velocities  mark  the  onset  of  the  nonlinear  regime,  which  is  dominated  by 
cusps  formed  in  the  tightly  packed  turbulent  flame.  The  nonlinear  increase  of  S  /  in  cusps,  and,  therefore,  the  nonlinear 
dependence  on  At  of  the  contribution  of  cusps  to  5  7-,  causes  a  strong  exaggerated  response  of  S  7-  to  the  increase  in 
At. 

The  onset  of  the  nonlinear  regime  of  flame  evolution  marks  the  breakdown  of  Damkohler’s  concept.  Moreover,  in 
this  regime,  flame  propagation  can  no  longer  be  viewed  as  a  local  process.  In  particular,  the  local  flame  speed  at  each 
point  of  the  flame  surface  is  no  longer  determined  only  by  the  local  thermodynamic  state  of  the  flow  or  by  turbulent 
motions  on  scales  A  <  6l .  Instead,  in  this  regime,  S  /  is  also  determined  by  long-range  velocity  correlations,  which 
produce  flame  collisions  and  can  span  the  full  size  of  the  system. 

Both  the  criteria  given  by  eqs.  (27)-(29)  and  the  results  of  numerical  simulations  show  that  the  transition  to  the 
nonlinear  regime  occurs  well  within  the  thin  reaction  zone  mode  of  combustion  and,  in  particular,  at  Ka  >  20  (Fig.  15). 
Above  this  critical  value  of  Ka,  the  flame  evolution  likely  remains  in  the  nonlinear  regime  at  all  turbulent  intensities, 
since  the  turbulence  is  more  efficient  at  packing  the  flame  than  at  broadening  it.  In  particular,  at  high  enough  values 
of  Ui ,  turbulence  will  eventually  break  the  internal  flame  structure,  which  will  increase  the  local  flame  speed.  This 
will,  in  turn,  make  the  right-hand  side  of  eq.  (27)  larger,  thereby,  enhancing  the  stabilizing  effect  of  cusps.  On  the 
other  hand,  the  flame  width.  If,  will  also  grow,  increasing  the  critical  flame  separation  necessary  for  the  onset  of  the 
nonlinear  regime.  Since  S  /  is  determined  by  smaller-scale  turbulence,  while  flame  folding  is  governed  by  the  faster 
larger-scale  motions,  it  is  unlikely  that  S /  can  grow  fast  enough  to  compensate  for  the  increase  in  both  If  and  the 
turbulent  speeds  which  fold  the  flame.  Consequently,  at  higher  U/,  i.e.,  even  in  the  broken  reaction  zone  mode,  the 
linear  regime  most  likely  cannot  be  recovered.  This  issue,  however,  requires  further  investigation  in  future  studies. 

The  increase  of  S  /  in  cusps  is  discussed  here  for  Le  —  1.  It  is  well  known,  however,  that  when  Le  4  1,5/  can 
increase  with  the  curvature  of  the  flame  [3],  Therefore,  due  to  the  large  flame  curvature  in  cusps,  any  imbalance  be¬ 
tween  thermal  and  diffusion  fluxes  can  significantly  exacerbate  the  enhancement  of  5  /.  Consequently,  the  exaggerated 
response  of  S  7-,  discussed  here,  can  be  substantially  larger  for  the  reactive  flows  characterized  by  Le  4  1. 

Results  presented  in  this  work  show  that  for  the  large  range  of  turbulent  intensities  and  system  sizes,  knowledge 
of  At  is  no  longer  sufficient  to  predict  the  magnitude  of  St-  For  instance,  at  the  turbulent  speeds  considered  here, 
using  At  as  a  guide  results  in  errors  on  the  order  of  30  -  50%.  Such  errors  are  quite  substantial  given  that  the  flow 
evolution  is  typically  very  sensitive  to  the  rate  of  energy  release.  Therefore,  it  is  particularly  important  to  account  for 
the  nonlinear  effects  in  subgrid-scale  models  that  primarily  focus  on  determining  the  evolution  of  At  (e.g.,  [7]). 

Furthermore,  the  formation  of  cusps  and  the  resulting  rapid  flame  propagation  in  certain  regions  is  a  crucial  part 
of  the  turbulent  flame-brush  evolution  in  the  high-speed  regime.  Thus,  properly  capturing  this  regime  in  numerical 
models  requires  the  domain  size  to  be  larger  than  the  integral  scale  in  order  to  accommodate  the  folding  of  the  flame 
by  turbulence.  Making  the  domain  smaller  than  /  would  significantly  hamper  this  process,  while  making  the  domain 
smaller  than  ~6l  would  completely  eliminate  it. 

Finally,  two  possibilities  exist  for  the  flame  evolution  in  the  nonlinear  regime.  If  the  turbulent  intensity  is  not 
too  high,  an  equilibrium  is  established  between  flame-surface  creation  and  its  rapid  destruction  in  cusps.  This  results 
in  the  turbulent  flame  propagating  in  a  steady  state,  which  is  the  situation  observed  here.  On  the  other  hand,  if  the 
turbulent  intensity  continues  to  rise,  eventually  Xy  and  At  will  approach  their  limiting  values  given  by  eqs.  (33)  -(34). 
In  this  regime,  as  was  discussed  in  §  5.6,  both  5  r  and  dSi/cIA/  can  become  arbitrarily  large.  Such  singular  behavior 
suggests  that,  in  reality,  at  some  point  the  steady  state  ceases  to  exist  and  the  system  must  undergo  a  qualitative 
transformation  in  order  to  accommodate  the  rise  in  S  t  or,  equivalently,  in  the  rate  of  energy  release  per  unit  volume. 
Such  qualitative  change  may  indicate  the  transition  from  the  deflagration  to  a  detonation.  Detailed  discussion  of  this 
non-steady  regime  will  be  presented  in  a  separate  paper. 
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Appendix  A.  Definition  of  the  flame  surface  density 


The  definition  of  Z(T),  given  in  eq.  (18),  relied  on  the  isosurface  area  normalization  by  the  full  width  of  the  flame 
brush.  It  is,  however,  not  clear  a  priori  whether  all  isosurfaces  occupy  the  full  interior  of  the  flame  brush  and  are  not 
confined  to  a  smaller  part  of  it.  Consider  Fig.  2.  It  can  be  seen  that  neither  the  Y  -  0.05,  nor  the  Y  =  0.95  isosurface 
extend  over  the  full  width  of  the  flame  brush,  and  they  indeed  occupy  a  smaller  volume. 

In  order  to  give  this  a  more  quantitative  representation,  we  define  ZQjnax  and  z\min  by  analogy  with  zo,mi„  and  z\max 
that  were  used  to  specify  5r  in  eqs.  (9)-(  10).  In  particular, 


—Ojucix  =  min(z)  :  Y(x,y,z)  >  0.05  V  (. x,y,z  >  z0,max), 
Zijnin  =  max(z) :  Y(x,y,z)  <  0.95  V  (x,  y,  z  <  zhmm). 


(35) 


In  other  words,  ZQ^max  is  the  z-coordinate  of  the  rightmost  cell  with  pure  product,  while  z\mm,  respectively,  is  the 
z-coordinate  of  the  leftmost  cell  with  pure  fuel.  Thereby,  Z(]x,,ax  and  z\  mn  effectively  measure  the  furthest  extent  of 
product  and  fuel  penetration  into  the  flame  brush.  Given  the  definitions  of  zojm;„  and  z\jnax,  provided  in  §  2.2,  zo,m,„ 
and  z§max  bound  the  volume  confining  the  Y  =  0.05  isosurface,  while  z\mn  and  z\  max  bound  the  Y  =  0.95  isosurface. 
At  the  same  time,  if  zi,„„„  <  ZQjnax,  these  two  variables  bound  the  region  of  macroscopic  mixing  of  product  and  fuel  in 
the  flame  brush,  i.e.,  the  region  in  which  both  pure  fuel  and  pure  product  can  be  found.  Fig.  2  provides  the  illustration 
of  all  four  of  these  quantities. 

Next,  we  define  the  position  of  the  turbulent  flame  brush  as 


zT 


(36) 


Modified  zo,maX  and  z\  min  are  then  defined  as  offsets  of  ZQmax  and  zijTO„  with  respect  to  zj-  normalized  by  6r,  namely 


-  _  "O.Mfl.V  ~  ZT  _  _  Z 1  min  -  ZT 

z0,max  —  ,  Z  \  min  —  ..  { Z>  / ) 

or  or 

Thus,  z§  max  and  zi,„„„  are  the  relative  measure  of  the  extent  to  which  pure  fuel  and  pure  product  penetrate  into  the 
flame  brush.  It  follows  from  the  definition  (37),  that  both  quantities  take  on  values  in  the  interval  [-0.5, 0.5].  For 
example,  in  a  planar  laminar  flame,  they  are  constant  with  z o,max  =  -0.5,  ziim;„  =  0.5.  In  the  turbulent  flame  brush, 
zero  values  would  indicate  that  both  fuel  and  product  reach  the  midpoint  of  the  flame  brush,  they  are  confined  to 
the  left  and  right  halves  of  the  brush,  and  there  is  no  macroscopic  mixing  of  fuel  and  product.  If  z^max  =  0.5  and 
zi  min  =  -0.5,  then  both  pure  product  and  pure  fuel  can  be  found  throughout  the  entire  volume  of  the  flame  brush. 

Fig.  16  shows  the  evolution  of  Z(imax  and  z\  mm  for  simulations  SI  -  S3.  In  all  cases,  both  parameters  exhibit  fairly 
similar  evolution  oscillating  around  zero.  Overall,  pure  fuel  and  product  tend  to  be  separated  in  the  flame  brush.  As 
the  system  evolves,  however,  it  undergoes  recurring  transitions  between  periods  of  enhanced  fuel-product  mixing  and 
episodes  of  their  near  complete  separation.  The  correlation  between  these  quantities  and  S  r  is  much  less  prominent 
compared,  for  instance,  with  that  between  A(Y)  and  Sr,  although  it  is  possible  to  associate  some  peaks  and  troughs 
with  the  corresponding  changes  in  the  turbulent  flame  speed. 

It  follows  from  Fig.  16  that  throughout  the  course  of  the  simulation,  both  the  Y  =  0.05  and  Y  =  0.95  isosurfaces 
occupy  the  region  smaller  than  the  full  flame-brush  volume.  This  is  also  the  case  for  other  isosurfaces  (cf.  Fig.  2). 
Consequently,  an  argument  can  be  made  that  when  calculating  Z.(Y ),  normalization  should  be  performed  not  over  the 
full  volume  of  the  flame  brush,  but  rather  over  the  respective  volume  bounding  the  given  isosurface.  For  instance,  for 
Y  =  0.05  this  volume  is  Lr(zojnax  -  ZQmln).  Such  definition  of  X(T)  would  be  a  more  accurate  measure  of  how  tightly 
a  given  isosurface  is  folded.  Fig.  16  shows  that,  on  average,  both  the  Y  =  0.05  and  Y  =  0.95  isosurfaces  are  confined 
to  about  half  of  the  total  volume  of  the  flame  brush.  Then  with  this  modified  definition  of  £,  its  corresponding  values 
in  Figs.  7-8  should  be  multiplied  by  a  factor  of  two.  Ultimately,  however,  we  are  interested  in  determining  how  the 
overall  flamelet,  rather  than  individual  isosurfaces,  is  folded  inside  the  flame  brush  since  only  the  flamelet  as  a  whole 
has  actual  physical  significance.  Therefore,  we  find  that  the  uniform  normalization  over  the  total  volume  occupied  by 
the  flamelet  represents  a  more  physically  grounded  choice. 

Further  motivation  for  our  choice  of  normalization  in  the  definition  of  E (18)  comes  when  one  considers  how  the 
values  of  zo,max  and  z\mn  would  change  for  larger  system  sizes.  Pure  fuel  and  product  are  always  separated  by  the  full 
width  of  the  flamelet.  In  particular,  in  our  system  the  characteristic  distance  between  the  points  Y  =  0.05  and  Y  =  0.95 
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in  the  flamelet  structure  is  &26l  (see  Fig.  7  in  [9]).  It  then  follows  that 


,max  z0 ,max  ^  2hy,  Z\  min  -O.niin  ^7  2S [  .  (38) 

Using  eqs.  (8),  (36),  and  (37),  inequalities  (38)  can  be  transformed  into  the  following  conditions  which  must  be 
satisfied  at  each  moment  in  time 


ZO  ,max 


<0.5  -¥■ 

Oj 


Zij 


>-0.5  +  ^. 

Oj 


(39) 


Horizontal  dashed  lines  in  Fig.  16  show  the  average  limiting  z$max  and  z\  mm  based  on  the  values  of  6t/Sl  listed  in 
Table  3  for  each  simulation.  In  particular,  in  all  three  calculations  zo,max  must  be  <  0.36  while  z\  min  must  be  >  -0.36, 
which  closely  agrees  with  the  data  shown  in  Fig.  16. 

This  demonstrates  that  both  isosurfaces  of  Y  —  0.05  and  Y  =  0.95  can  never  occupy  the  full  volume  of  the  flame 
brush  simply  because,  by  definition,  they  are  separated  from  the  boundary  of  the  flame  brush  by  the  flamelet  thickness. 
The  latter,  in  our  case,  is  a  substantial  fraction  of  the  total  flame  brush  width.  It  was  discussed  in  §  5.6  that  increases 
with  the  turbulent  integral  scale,  or,  equivalently,  with  the  system  size  and  energy  injection  scale  ([9],  also  see  [10]). 
It  then  follows  from  eq.  (39)  that  as  5t/6l  — >  oo,  then  zotfnax  — >  0.5  and  z\  mln  — »  -0.5.  Therefore,  in  larger  systems  the 
flame  brush  width  becomes  large  in  comparison  with  the  width  of  the  flamelet  and,  as  a  result,  the  volume  bounding 
each  isosurface  becomes  well  approximated  by  the  total  volume  of  the  flame  brush.  Consequently,  in  the  limit  of  large 
values  of  6t,  the  uniform  normalization  for  all  values  of  Y  in  eq.  (18)  becomes  equivalent  to  the  normalization  by  the 
actual  volume  bounding  a  given  isosurface. 
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Table  1 .  Input  model  parameters  and  resulting  computed  laminar  flame  properties 


Input 


To 

293  K 

Po 

1.01  x  106  erg/cm3 

Po 

8.73  x  10-4  g/cm3 

7 

1.17 

M 

21  g/rnol 

B 

6.85  x  1012  cm3/(g  s) 

Q 

46.37  RT0 

<7 

43.28  RT0  /  M 

Ko 

2.9  x  10~5  g/(s  cm  Kn) 

Do 

2.9  x  10~5  g/(s  cm  Kn) 

n 

0.7 

Output 

TP 

2135  K 

Pp 

1.2  x  10~4  g/cm3 

6l 

0.032  cm 

Sl 

302  cm/s 

Initial  temperature 
Initial  pressure 
Initial  density 
Adiabatic  index 
Molecular  weight 
Pre-exponential  factor 
Activation  energy 
Chemical  energy  release 
Thermal  conduction  coefficient 
Molecular  diffusion  coefficient 
Temperature  exponent 

Post-flame  temperature 
Post-flame  density 
Laminar  flame  thermal  width 
Laminar  flame  speed 


Table  2.  Parameters  of  simulations3 


SI 

S2 

S3 

Description 

£> 

64  x  64  x  1024 

128  x  128  x  2048 

256  x  256  x  4096 

Domain  grid  size 

Da 

1  x  1  x  16 

Domain  aspect  ratio 

L 

0.259  cm  =  8 SL 

Domain  width,  energy-injection  scale 

A.v 

4.05  x  10~3  cm 

2.02  x  10~3  cm 

1.01  x  10~3  cm 

Cell  size 

8 

16 

32 

<5l/A.v 

ZT,0 

1.95  cm  =  7.52L 

Initial  flame  position  along  z-axis 

£ 

1.26  x  109  erg/(cm3  s) 

Energy-injection  rate 

Us 

4.53  x  103  cm/s  =  15SL 

Turbulent  velocity  at  scale  Sl 

u 

9.07  x  103  cm/s  =  305  L 

Turbulent  velocity  at  scale  L 

U rms 

1.04  x  104  cm/s  =  34.485  L 

Turbulent  r.m.s.  velocity 

u, 

5.60  x  103  cm/s  =  18.545 L 

Integral  velocity 

i 

6.04  x  10~2  cm  =  1 .8767 

Integral  scale 

Ted 

2.86  x  10~5  s 

Eddy  turnover  time,  L/U 

hgn 

3.0red 

3.0  Ted 

2.0red 

Time  of  ignition 

t total 

16.0  Ted 

Total  simulation  time  from  tign 

Da 

0.05 

Damkohler  number?  (1/U/)/(2Sl/S l) 

La 

9.47  x  10~6  cm  =  2.96  x  10~4<5z. 

Gibson  scale,  l(S  l/Ui)3 

Map 

0.25 

Mach  number  in  fuel,  U(-yPo/po)  1/2 

MaP 

0.09 

Mach  number  in  product,  U(yPo!pp) ~1/2 

a  Parameters  common  to  all  simulations  are  shown  only  once  in  S2  column. 
b  In  this  definition  of  Da,  2 Sl  indicates  the  full  flame  width.  See  [9]  for  further  discussion. 
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Table  3.  Time-averaged  properties  of  the  turbulent  flame  brush3 


St/6l 

0(ST/SL) 

S  t/S  L 

0(ST/SL) 

Aq.is/L2 

O(A0A5/L2) 

S015,  mm  1 

CXSo.is) 

fo.15  000.15) 

SI 

16.13 

6.09 

3.91 

1.52 

S2 

14.86 

1.96 

4.50 

2.29 

3.23 

2.84 

4.07 

1.39  1.29 

S3 

14.42 

4.09 

3.12 

1.30 

Time-averaging  for  all  variables  is  performed  over  the  time  interval  [2 red  -  16rej], 
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Figure  1  (a)  Evolution  of  the  turbulent  flame  width  6t  normalized  by  ()/,.  Note  that  domain  width  L  =  8 5l  indicated 
with  the  horizontal  dashed  line,  (b)  Evolution  of  the  turbulent  flame  speed  5  t  normalized  by  S  l.  In  both  panels:  black 
lines  correspond  to  simulation  SI,  red  to  S2,  and  green  to  S3.  (Reproduced  from  [9].) 


Figure  2  Isosurfaces  of  Y  in  simulation  S2  at  t  =  1 3 Tetj  (cf.  Fig. 3,  middle  row,  left  panel  in  [9]).  Isosurface  values 
are  0.05  (red),  0.6  (green),  0.95  (blue).  Red  and  green  isosurfaces  bound  the  flamelet  reaction  zone.  Green  and  blue 
isosurfaces  bound  the  preheat  zone.  The  zo,m;„  and  z\^max  mark  the  flame-brush  bounds.  The  z^max  and  z\mi„  indicate, 
respectively,  the  maximum  extents  of  product  and  fuel  penetration  into  the  flame  brush  (see  Appendix  A  for  further 
discussion  of  Z(imax  and  Z|., ,„■„).  (Reproduced  from  [9].) 
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Figure  3  Evolution  of  the  normalized  flame  surface  area  (a)  zlo.oi/E2  and  (b)  A0.99/L2.  (c)  Time-averaged  normalized 
flame  surface  area  A(Y)/L2.  Time  averaging  is  performed  over  the  time  interval  [2 red  -  1 6r«/],  circles  represent 
calculated  values  and  solid  lines  are  the  Akima  spline  fits.  Shaded  gray  region  shows  the  distribution  of  the  reaction 
rate,  Y,  in  the  exact  laminar  flame  solution  nonnalized  by  its  peak  value  Ymax  =  9.5  x  104  s  1 .  In  all  panels:  black 
corresponds  to  simulation  SI,  red  to  S2,  and  green  to  S3. 


Figure  4  (a)  Evolution  of  the  stretch  factor  /0.15.  (b)  Time-averaged  distribution  of  the  stretch  factor  I(  Y ).  Circles 
represent  calculated  values  and  solid  lines  are  the  Akima  spline  fits.  Time  averaging  is  performed  over  the  time 
interval  [2 red  -  1  6t,,,/].  The  shaded  gray  area  shows  +<x  standard  deviation  of  the  instantaneous  values  of  l(  Y )  in 
simulation  S3.  In  both  panels:  black  corresponds  to  simulation  SI,  red  to  S2,  and  green  to  S3. 


26 


7.0 


Figure  5  Correlation  between  A(Y)/L 2  and  S -/  /S /  in  simulation  S3.  A(Y)  is  considered  for  three  values  of  Y :  (a,  b) 
Y  =  0.15,  (c,  d)  Y  -  0.5,  and  (e,  f)  Y  =  0.95.  In  panels  (a),  (c),  and  (e)  A/L2  is  shifted  in  time  to  the  right  by  the  time 
lag  A tc  shown  in  each  panel.  Graphs  in  panels  (b),  (d),  and  (f)  are  constructed  based  on  the  data  shown,  respectively, 
in  panels  (a),  (c),  and  (e)  excluding  the  time  (0  -  2 )Ted  during  which  the  turbulent  flame  develops  its  equilibrium  state. 
In  panel  (b)  each  data  point  is  colored  according  to  the  corresponding  value  of  the  surface  density  L(T)  with  the  color 
scale  given  in  units  mm  1  (cf.  Fig.  8a, b).  In  panels  (b),  (d),  and  (f)  dashed  lines  show  time-averaged  values  of  A/L2 
and  St/S l,  solid  lines  show  the  least  squares  fit  with  its  slope  given  in  each  panel,  and  dash-dot  lines  correspond  to 
St/S  l  =  A/L2.  Shaded  gray  regions  in  panels  (a)-(e)  illustrate  the  exaggerated  response  of  Sr  to  the  increase  of  A. 
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Figure  6  Correlation  between  A(Y)/L 2  and  6j/L  in  simulation  S3.  A(Y )  is  considered  for  three  values  of  V :  (a,  b) 
Y  =  0.15,  (c,  d)  Y  =  0.5,  and  (e,  f)  Y  =  0.95.  Graphs  in  panels  (b),  (d),  and  (f)  are  constructed  based  on  the  data 
shown,  respectively,  in  panels  (a),  (c),  and  (e)  excluding  the  time  (0  -  2)r„/  during  which  the  turbulent  flame  develops 
its  equilibrium  state.  In  panels  (b),  (d),  and  (f)  dashed  lines  show  time-averaged  values  of  A/L 2  and  5j/L,  while  the 
solid  lines  show  the  least  squares  fit  with  its  slope  given  in  each  panel. 
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Figure  7  Evolution  of  the  flame  surface  density  (a)  Xo.oi  and  (b)  X0.99.  (c)  Time-averaged  flame  surface  density  X(T). 
Time  averaging  is  performed  over  the  time  interval  [2 Tej  -  1 6rfvy],  circles  represent  calculated  values  and  solid  lines 
are  the  Akima  spline  fits.  Shaded  gray  region  shows  the  distribution  of  the  reaction  rate,  Y,  in  the  exact  laminar  flame 
solution  normalized  by  its  peak  value  Ymax  =  9.5  x  104  .  In  all  panels:  black  corresponds  to  simulation  SI,  red  to 

S2,  and  green  to  S3. 
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Figure  8  Correlation  between  X(  Y )  and  St/Sl  in  simulation  S3.  1(b)  is  considered  for  three  values  of  F:  (a,  b) 
Y  =  0.15,  (c,  d)  F  =  0.5,  and  (e,  f)  Y  =  0.95.  In  panels  (a),  (c),  and  (e)  X  is  shifted  in  time  to  the  right  by  the  time  lag 
Ac  shown  in  each  panel.  Graphs  in  panels  (b),  (d),  and  (f)  are  constructed  based  on  the  data  shown,  respectively,  in 
panels  (a),  (c),  and  (e)  excluding  the  time  (0  -  2)t(,j  during  which  the  turbulent  flame  develops  its  equilibrium  state. 
In  panels  (b),  (d)  and  (f)  dashed  lines  show  time-averaged  values  of  X  and  St/Sl ,  while  the  solid  lines  show  the  least 
squares  fit  with  its  slope  given  in  each  panel. 
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Figure  9  Flame  collision  and  the  formation  of  cusps  in  the  turbulent  flame.  Shown  is  the  flame-brush  structure  based 
on  the  isovolume  of  the  fuel  mass  fraction  in  simulation  S3.  Bounding  isosurfaces  represent  Y  =  0.05  and  Y  =  0.95 
and  the  flame  brush  is  shown  from  the  product  side.  Upper  panel  corresponds  to  the  time  t  =  I  1 .86red,  while  the 
middle  and  lower  panels  show  the  flame  structure,  respectively,  0.1  Ted  =  3  /is  and  0.2t(Y/  =  6  /us  later.  The  thin  black 
line,  corresponding  to  Y  —  0.6,  marks  the  boundary  between  the  preheat  and  reaction  zones.  The  thin  white  line, 
corresponding  to  Y  —  0.2,  shows  the  location  of  the  peak  reaction  rate.  Regions  A,  B,  and  C  show  the  three  main 
stages  of  the  flame  collision  and  the  formation  of  a  cusp  discussed  in  §  5.2.  Note  also  two  elongated  regions  of  flame 
collision  forming  near  the  upper  face  of  the  domain. 
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Figure  10  Structure  of  a  cusp  formed  by  the  collision  of  parallel  flame  sheets.  Shown  is  the  distribution  of  Y  for  the 
flame  inclination  angle  a  =  1°  (upper  panel)  and  a  =  4°  (lower  panel).  Scale  of  the  panel  axes  is  given  in  units  of 
5l-  Away  from  the  cusp  tip,  the  flame  propagates  in  the  direction  normal  to  its  surface  with  the  laminar  flame  speed, 
S l,  causing  the  tip  to  move  to  the  right  with  the  speed  Dc.  Also  shown  is  the  length,  lc,  of  the  collision  region  (see 
text).  Thin  black  line  marks  the  boundary  between  the  reaction  and  preheat  zones,  while  the  thin  white  line  indicates 
the  region  of  peak  reaction  rate.  Note  the  substantial  broadening  of  the  flame,  and  thus  the  reaction  zone,  near  the  tip 
of  the  cusp  at  a  —  1°  compared  to  a  =  4°  (cf.  Fig.  9). 
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Figure  1 1  Structure  of  a  cusp  formed  by  the  collision  of  parallel  flame  sheets  for  four  flame  inclination  angles,  a. 
Shown  are  distributions  of  (a)  the  fuel  mass  fraction,  Y,  and  temperature,  T,  as  well  as  (b)  the  reaction  rate,  Y ,  along 
the  symmetry  line  of  the  cusp  (cf.  Fig.  10).  Dashed  lines  in  both  panels  indicate  the  exact  planar  laminar  flame 
solution.  Profiles  of  T  and  Y  are  normalized  by  their  respective  maximum  values  in  the  exact  laminar  solution,  i.e., 
TP  (see  Table  1 )  and  Ymax  =  9.5  x  104  s'1. 
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Figure  12  (a)  Dependence  of  the  normalized  speed  of  cusp  propagation,  Dc/Sl,  on  the  flame  inclination  angle,  a. 
Solid  line  corresponds  to  the  analytic  expression  given  in  eq.  (20),  red  circles  show  the  computed  values.  Horizontal 
dashed  line  indicates  sound  speed  in  cold  fuel,  (b)  Dependence  of  the  maximum  normalized  local  burning  speed  in 
the  cusp,  S  * /S  given  by  eq.  (21),  on  a.  Red  circles  show  the  computed  values,  solid  line  is  the  Akima  spline  fit. 


Figure  13  Dependence  of  the  stretch  factor,  I,  given  by  eq.  (23)  on  the  nonnalized  length  of  the  flame-collision  region, 
IJ6l,  for  four  values  of  a.  Shaded  gray  area  indicates  the  full  laminar  flame  width  2c')/  .  See  text  for  further  details. 


Figure  14  Illustration  of  the  idealized  perturbed  flame  stabilized  by  the  propagation  of  a  cusp  with  the  speed  Dc  (cf. 
Fig.  10). 
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Figure  15  Combustion  regime  diagram  according  to  [10].  Orange  region  above  the  Ka  ~  20  line  shows  the  range  of 
the  regimes  in  which  the  formation  of  cusps  is  expected  to  result  in  values  of  I  substantially  above  unity  according 
to  the  criterion  discussed  in  §  5.5.  Filled  red  square  corresponds  to  the  simulations  presented  in  this  work,  while  the 
open  red  square  shows  the  regime  in  which  the  value  of  Z0.15  ~  1.16  was  determined.  Flamelets  are  typically  suggested 
to  exist  in  the  regimes  below  the  K as  =  1  line  [10].  The  traditional  form  of  the  diagram  was  also  modified  by  adding 
the  Map  —  U(yPo/po)  1,2  =  1  line  indicating  the  region  of  supersonic  turbulence  in  the  cold  H2-air  fuel  under  the 
atmospheric  conditions.  See  text  for  further  details. 
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Figure  16  Evolution  of  the  normalized  maximum  extents  of  product  and  fuel  penetration  into  the  flame  brush  Z()  max 
and  z |  min  (see  text  for  the  definition  and  Fig.  2  for  the  illustration)  for  simulation  SI  (a),  S2  (b),  and  S3  (c).  ZQmax  is 
shown  with  red  line,  z\mln  -  with  blue  line.  Shaded  regions  mark  the  extent  of  macroscopic  mixing  of  pure  fuel  and 
product  inside  the  flame  brush.  Horizontal  dashed  lines  show  the  average  limiting  values  of  ZQmax  and  z\  min  given  by 
the  condition  (39)  and  calculated  using  the  values  of  <)■/  / Sl  listed  in  Table  3. 
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